The first two chunks of this r markdown file after the r setup allow for plot zooming, but it also means that the html file must be opened in a browser to view the document properly. When it knits in RStudio the preview will appear empty but the html when opened in a browser will have all the info and you can click on each plot to Zoom in on it.

Before you begin

Notes

A few notes about this script.

If you are running this with the 2022-2023 data make sure you download the whole (OSM_2022-2023 GitHub repository)[https://github.com/ACMElabUvic/OSM_2022-2023] from the ACMElabUvic GitHub. This will ensure you have all the files, data, and proper folder structure you will need to run this code and associated analyses.

Also make sure you open RStudio through the R project (OSM_2022-2023.Rproj) this will automatically set your working directory to the correct place (wherever you saved the repository) and ensure you don’t have to change the file paths for some of the data.

Lastly, if you are looking to adapt this code for a future year of data, you will want to ensure you have run the ACME_camera_script_9-2-2024.R or .Rmd with your data as there is much data formatting, cleaning, and restructuring that has to be done before this code will work.

If you have question please email the most recent author, currently

Marissa A. Dyck
Postdoctoral research fellow
University of Victoria
School of Environmental Studies
Email: marissadyck17@gmail.com

Install packages

If you don’t already have the following packages installed, use the code below to install them.

install.packages('tidyverse') 
install.packages('ggpubr')
install.packages('corrplot')
install.packages('Hmisc')
install.packages('glmmTMB')
install.packages('MuMIn')
install.packages('TMB', type = 'source')
install.packages('rphylopic')

Load libraries

Then load the packages to your library.

library(tidyverse) # data tidying, visualization, and much more; this will load all tidyverse packages, can see complete list using tidyverse_packages()
library(ggpubr) # make modificaions to plot for publication (arrange plots)
library(PerformanceAnalytics)    #Used to generate a correlation plot
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## ################################### WARNING ###################################
## # We noticed you have dplyr installed. The dplyr lag() function breaks how    #
## # base R's lag() function is supposed to work, which breaks lag(my_xts).      #
## #                                                                             #
## # Calls to lag(my_xts) that you enter or source() into this session won't     #
## # work correctly.                                                             #
## #                                                                             #
## # All package code is unaffected because it is protected by the R namespace   #
## # mechanism.                                                                  #
## #                                                                             #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## # You can use stats::lag() to make sure you're not using dplyr::lag(), or you #
## # can add conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop   #
## # dplyr from breaking base R's lag() function.                                #
## ################################### WARNING ###################################
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
library(Hmisc) # used to generate histograms for all variables in data frame
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(glmmTMB)      #Constructing GLMMs
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.4.1
## Current Matrix version is 1.5.3
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
library(MuMIn) # for model selection
library(rphylopic) # add animal silhouettes to graphs
## You are using rphylopic v.1.3.0. Please remember to credit PhyloPic contributors (hint: `get_attribution()`) and cite rphylopic in your work (hint: `citation("rphylopic")`).

Data

Load detection data

Read in saved and cleaned detection data from the ACME_camera_script_9-2-2024.R.

# detection data
# read in saved and cleaned detection data from the ACME_camera_script_9-2-2024.R 
detections <- read_csv('data/processed/OSM_2022_ind_det.csv') %>% 
  
  # change site, species and event_id to factor
  mutate_if(is.character,
            as.factor)
## Rows: 14102 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (4): array, site, species, event_id
## dbl  (3): month, year, timediff
## dttm (1): datetime
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Data formatting

In order to get plots that have the same formatting as last years’ report we have to do a bit of data formatting. First we need to make sure we are including the same relevant species (some were ignored for last years’ report or grouped together)

Last years report had the following species

  • white-tailed deer
  • snowshoe hare
  • black bear
  • coyote
  • red squirrel
  • fisher
  • unknown
  • moose
  • lynx
  • spruce grouse
  • red fox
  • striped skunk
  • ruffed grouse
  • owl
  • grey wolf
  • domestic dog
  • cougar
  • raven
  • other
  • mule deer

And they grouped all humans except for staff as ‘Humans’. Let’s look at the species we have in this year’s data and try to format it the same way

detections %>% 
  
  # group by array and species
  group_by(array, species) %>% 
  summarise(n = n()) %>% 
  
  # have R print everything
  print(n = nrow(.))
## `summarise()` has grouped output by 'array'. You can override using the
## `.groups` argument.
## # A tibble: 119 × 3
## # Groups:   array [4]
##     array species                 n
##     <fct> <fct>               <int>
##   1 LU01  Beaver                  1
##   2 LU01  Black bear            380
##   3 LU01  Cougar                  7
##   4 LU01  Coyote                581
##   5 LU01  Domestic dog            6
##   6 LU01  Fisher                111
##   7 LU01  Grey jay               14
##   8 LU01  Grey wolf              21
##   9 LU01  Human                   3
##  10 LU01  Lynx                   55
##  11 LU01  Moose                  99
##  12 LU01  Other                   1
##  13 LU01  Other birds            60
##  14 LU01  Otter                   2
##  15 LU01  Owl                     2
##  16 LU01  Porcupine               5
##  17 LU01  Raven                   6
##  18 LU01  Red fox                50
##  19 LU01  Red squirrel          879
##  20 LU01  Ruffed grouse          14
##  21 LU01  Short-tailed weasel     5
##  22 LU01  Snowshoe hare        1443
##  23 LU01  Spruce grouse          12
##  24 LU01  Staff                  71
##  25 LU01  Striped skunk          39
##  26 LU01  Unknown               210
##  27 LU01  Unknown canid          48
##  28 LU01  Unknown deer          175
##  29 LU01  Unknown mustelid       13
##  30 LU01  Unknown ungulate        8
##  31 LU01  White-tailed deer    1953
##  32 LU13  ATVer                  31
##  33 LU13  Black bear            275
##  34 LU13  Caribou                 3
##  35 LU13  Coyote                187
##  36 LU13  Fisher                  5
##  37 LU13  Grey jay                2
##  38 LU13  Grey wolf              52
##  39 LU13  Human                   2
##  40 LU13  Hunter                  1
##  41 LU13  Long-tailed weasel      1
##  42 LU13  Lynx                  115
##  43 LU13  Marten                 27
##  44 LU13  Moose                 128
##  45 LU13  Other birds            12
##  46 LU13  Owl                     1
##  47 LU13  Red fox                 2
##  48 LU13  Red squirrel          240
##  49 LU13  Ruffed grouse           7
##  50 LU13  Short-tailed weasel     7
##  51 LU13  Snowshoe hare         573
##  52 LU13  Spruce grouse          25
##  53 LU13  Staff                  82
##  54 LU13  Striped skunk           1
##  55 LU13  Unknown                86
##  56 LU13  Unknown canid          10
##  57 LU13  Unknown deer            5
##  58 LU13  Unknown mustelid        3
##  59 LU13  White-tailed deer      86
##  60 LU13  Wolverine               8
##  61 LU15  ATVer                   1
##  62 LU15  Beaver                  2
##  63 LU15  Black bear            220
##  64 LU15  Canada goose            3
##  65 LU15  Caribou                51
##  66 LU15  Coyote                171
##  67 LU15  Fisher                 25
##  68 LU15  Grey jay               21
##  69 LU15  Grey wolf              61
##  70 LU15  Long-tailed weasel     15
##  71 LU15  Lynx                  122
##  72 LU15  Marten                 63
##  73 LU15  Moose                 157
##  74 LU15  Other birds            59
##  75 LU15  Otter                   5
##  76 LU15  Owl                     1
##  77 LU15  Red fox                39
##  78 LU15  Red squirrel          643
##  79 LU15  Ruffed grouse          11
##  80 LU15  Short-tailed weasel     7
##  81 LU15  Snowmobiler             1
##  82 LU15  Snowshoe hare         611
##  83 LU15  Spruce grouse          21
##  84 LU15  Staff                  78
##  85 LU15  Unknown                98
##  86 LU15  Unknown canid           7
##  87 LU15  Unknown deer           47
##  88 LU15  Unknown mustelid       16
##  89 LU15  Unknown ungulate        5
##  90 LU15  White-tailed deer     429
##  91 LU21  Black bear            544
##  92 LU21  Canada goose            1
##  93 LU21  Caribou                16
##  94 LU21  Cougar                  2
##  95 LU21  Coyote                 51
##  96 LU21  Fisher                 46
##  97 LU21  Grey jay               13
##  98 LU21  Grey wolf              55
##  99 LU21  Long-tailed weasel      1
## 100 LU21  Lynx                   72
## 101 LU21  Marten                 50
## 102 LU21  Moose                 233
## 103 LU21  Other                   1
## 104 LU21  Other birds            44
## 105 LU21  Owl                     8
## 106 LU21  Red fox                14
## 107 LU21  Red squirrel          219
## 108 LU21  Ruffed grouse          11
## 109 LU21  Short-tailed weasel     2
## 110 LU21  Snowmobiler             6
## 111 LU21  Snowshoe hare         284
## 112 LU21  Spruce grouse          19
## 113 LU21  Staff                  71
## 114 LU21  Unknown               162
## 115 LU21  Unknown canid           5
## 116 LU21  Unknown deer           65
## 117 LU21  Unknown mustelid       23
## 118 LU21  Unknown ungulate        4
## 119 LU21  White-tailed deer     839
# now let's create a new data frame (tibble) to work with for the OSM figure summaries specifically

# I personally would lump all the unknown together and all the birds together but for the sake of consistency with last years' figures we will remove some entries, let's create a vector of entries to drop

species_drop <- c('Staff',
                  'Unknown deer',
                  'Unknown ungulate',
                  'Unknown canid',
                  'Unknown mustelid',
                  'Other birds')

# now we can create the new data frame with some changes consistent w/ choices made for 2021-2022
detections <- detections %>% 
  
  # for summarizing, lets lump all the recreational humans into "Humans"
  mutate(species = recode_factor(species, 
                                 "Snowmobiler" = "Human",
                                 "ATVer" = "Human",
                                 'Hunter' = 'Human')) %>% 
  
  # remove species we don't want to plot
  filter(!species %in% species_drop)

We will also want to subset the data by landscape unit (LU) and generate a new data frame for each LU to use for plotting

# we will also want to create a data frame for each LU to plot individually

# LU1
dets_LU1 <- detections %>% 
  filter(array == 'LU01')

# LU13
dets_LU13 <- detections %>% 
  filter(array == 'LU13')

# LU15
dets_LU15 <- detections %>% 
  filter(array == 'LU15')

# LU21
dets_LU21 <- detections %>% 
  filter(array == 'LU21')

ANDREW HELP

Can you make the above code into a forloop which assigns each new data frame created from subsetting as dets_LUname?

Detection plots

Detection data

Now we can apply the same data formatting for each LUs’ data frame using purrr.

We want to count the number of independent detections per species per LU to use in the detection plots

# apply the same formatting to each LU data frame using purrr map
detection_data <- list(dets_LU1,
                       dets_LU13,
                       dets_LU15,
                       dets_LU21) %>% 
  
  purrr::map(
    ~.x %>% 
      
      # group by species
      group_by(species) %>%
      
      # calculate a column with unique accounts of each species
      mutate(count = n_distinct(event_id)) %>% 
      
      # keep just the columns we need
      select(species, count) %>% 
      
      # keep only unique (distinct) rows so we should be left with one row per species, this helps with plotting later if you don't do it ggplot will try to count and plot each row it's annoying
      distinct()) %>% 
  
  # set names of list objects
  purrr::set_names('Detections LU01',
                   'Detections LU13',
                   'Detections LU15',
                   'Detections LU21')

Detection plots

Now to graph independent detections for each LU using purrr, this avoids a TON of code repetition needed to plot each one individually

We use purrr::imap() instead of purrr::map() because imap maintains the variable names in our list (e.g. Detections LU01, Detections LU13, etc.) which we can then use to title each plot.

Within purrr::imap() we just paste the code we would use for a single ggplot since all the graphical elements (except the title which we change with the file name [.y]) are the same

# create object detection plots which uses the detection_data list (w/ all 4 LUs)
detection_plots <- detection_data %>% 
  
  # use imap instead of map as it allows us to use .y to paste the list element names as the plot titles later
  purrr::imap(
    ~.x %>% 
      
      # now just copy and paste the ggplot code for the detection graphs
      ggplot(.,
             aes(x = reorder(species, count), y = count)) +
      
      # plot as bar graph using geom_col so we don't have to provide a y aesthetic
      geom_col() +
      
      # switch the x and y axis
      coord_flip() +
      
      # add the number of detections at the end of each bar
      geom_text(aes(label = count),
                color = "black",
                size = 3,
                hjust = -0.3,
                vjust = 0.2) +
      
      # label x and y axis with informative titles
      labs(x = 'Species',
           y = 'Number of Independent (30 min) Detections') +
      
      # add title to plot with LU name the .y will take the name of whatever you named each list element in the detection_data list, so make sure this name is what you want on the ggtitle
      ggtitle(.y) +
      
      # set the theme
      theme_classic() +
      theme(plot.title = element_text(hjust = 0.5)))

# view plots, this will print each in it's own window so you have to scroll back in the plot viewer pane to look at each one
detection_plots
## $`Detections LU01`

## 
## $`Detections LU13`

## 
## $`Detections LU15`

## 
## $`Detections LU21`

Save detection plots

Now we want to save these plots in case we need each individual one (we will combine the detection and naive occ plots into a single figure for each LU later and use those for the OSM report, but we may want these standalone plots later so let’s save them while they are here).

We can save all the plots from the purrr iteration above using purrr::imap. imap is used instead of map because it allows us to retain the list object names (plot names) to paste as the file name with the .y command.

IMPORTANT if you are using this code for a future github repo, DO NOT use .tiff as the file extension. This will cause issues when trying to push any changes to the github repo as the files are too large to meet githubs requirements

# save plots only use if needed
purrr::imap(
  detection_plots,
  ~ggsave(.x,
             file = paste0("figures/",
                           .y,
                           '.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
          dpi = 600,
          width = 11,
          height = 9,
          units = 'in'))
## $`Detections LU01`
## [1] "figures/Detections LU01.jpg"
## 
## $`Detections LU13`
## [1] "figures/Detections LU13.jpg"
## 
## $`Detections LU15`
## [1] "figures/Detections LU15.jpg"
## 
## $`Detections LU21`
## [1] "figures/Detections LU21.jpg"

Naive occupancy

Data

We also need to alter the detection data a bit to use for naive occupancy plots.

We will use the individual LU detection data like we did before and use purrr::map() to apply the dame data formatting to all 4 data frames.

Here we want to calculate the total number of sites in each LU, the number of sites each species was detected at in each LU and then use both those numbers to calculate naive occupancy for each species in each LU

# First we need to alter the data frame a bit for these plots, let's create a data frame for each LU (I couldn't figure out how to do this without assigning individual data frames for each UGH)


# apply the same formatting to each data frame using purrr
occupancy_data <- list(dets_LU1,
                       dets_LU13,
                       dets_LU15,
                       dets_LU21) %>% 
  
  purrr::map(
    ~.x %>% 
      
      # calculate the total number of sites for each LU
      mutate(total_sites = n_distinct(site)) %>% 
      
      # group by species to calculate the number of sites each spp occurred at
      group_by(species) %>% 
  
      # add columns to count the number of sites each spp occurred at and then the naive occupancy
  reframe(count = n_distinct(site),
          naive_occ = count/total_sites,
          ind_det = n_distinct(event_id)) %>% 
  
    # keep just the columns we need
  select(species, naive_occ, ind_det) %>% 
  
    # keep only unique (distinct) rows so we should be left with one row per species, this helps with plotting
  distinct()) %>% 
  
  purrr::set_names('Naive Occupancy LU01',
                   'Naive Occupancy LU13',
                   'Naive Occupancy LU15',
                   'Naive Occupancy LU21')

Occupancy plots

Now we can graph naive occupancy for each LU using purrr, and as with the detection plots this saves a massive amount of coding using purrr to run an iteration on the data files and produce four plots at once instead of copying and pasting code for each individually

# create object occupancy_plots which uses the occupancy_data list (w/ all 4 LUs)
occupancy_plots <- occupancy_data %>% 
  
  # use imap instead of map as it allows us to use .y to paste the list element names as the plot titles later
  purrr::imap(
    ~.x %>% 

      # now just copy and paste the ggplot code for the occupancy graphs
      ggplot(.,
             aes(x = fct_reorder(species,
                                 ind_det), # this reorders the species so they match the order of the detection plot which makes it better for viewing when the plots are arranged together in 1 figure for each LU
                 y = naive_occ)) +
      
      # plot as bars using geom_col() which uses stat = 'identity', instead of geom_bar() which will count the rows in each group and plot that instead of naive occ
      geom_col() +
      
      # flip x and y axis 
      coord_flip() +
      
      # add text to end of bars that provides naive occ value
      geom_text(aes(label = round(naive_occ, 2)), 
                size = 3, 
                hjust = -0.3, 
                vjust = 0.2) +
      
      # relabel x and y axis and title
      labs(x = 'Species',
           y = 'Proportion of Sites With At Least One Detection') +
      
      # set plot title using .y (name of list object)
      ggtitle(.y) +
      
      # set. theme elements
      theme_classic()+
      theme(plot.title = element_text(hjust = 0.5)))

# view plots
occupancy_plots
## $`Naive Occupancy LU01`

## 
## $`Naive Occupancy LU13`

## 
## $`Naive Occupancy LU15`

## 
## $`Naive Occupancy LU21`

Save occupancy plots

As with the detection plots, we might want these individual plots later for something so we can use purrr::imap() to save them to the figures folder

Again avoid using the .tiff extension in github

# save plots 
purrr::imap(
  occupancy_plots,
  ~ggsave(.x,
          file = paste0("figures/",
                        .y,
                        '.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
          dpi = 600,
          width = 11,
          height = 9,
          units = 'in'))
## $`Naive Occupancy LU01`
## [1] "figures/Naive Occupancy LU01.jpg"
## 
## $`Naive Occupancy LU13`
## [1] "figures/Naive Occupancy LU13.jpg"
## 
## $`Naive Occupancy LU15`
## [1] "figures/Naive Occupancy LU15.jpg"
## 
## $`Naive Occupancy LU21`
## [1] "figures/Naive Occupancy LU21.jpg"

Final combined plots for report

The previous year’s report had a figure for each LU with the detections plot on the top and the occupancy plot on the bottom so we will recreate these for this year using ggarrange().

Unfortunately I could not figure out how to do this in purrr to reduce coding but luckily it isn’t too much repitition

# not sure I know how to do the following section in purrr just yet, but we've saved a ton of coding so far and it doesn't take much to arrange each of these individually

# LU1

# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU1_det_occ_plots <- ggarrange(detection_plots$`Detections LU01`, occupancy_plots$`Naive Occupancy LU01`,
                               labels = c("A", "B"),
                               nrow = 2)

# view plot
LU1_det_occ_plots

# LU13

# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU13_det_occ_plots <- ggarrange(detection_plots$`Detections LU13`, occupancy_plots$`Naive Occupancy LU13`,
                               labels = c("A", "B"),
                               nrow = 2)

# view plot
LU13_det_occ_plots

# LU15

# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU15_det_occ_plots <- ggarrange(detection_plots$`Detections LU15`, occupancy_plots$`Naive Occupancy LU15`,
                                labels = c("A", "B"),
                                nrow = 2)

# view plot
LU15_det_occ_plots

# LU21

# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU21_det_occ_plots <- ggarrange(detection_plots$`Detections LU21`, occupancy_plots$`Naive Occupancy LU21`,
                                labels = c("A", "B"),
                                nrow = 2)

# view plot
LU21_det_occ_plots

We can however, save all the figures again using purrr

# save all figures at once using purrr
final_det_occ_plots <- list(LU1_det_occ_plots,
                            LU13_det_occ_plots,
                            LU15_det_occ_plots,
                            LU21_det_occ_plots) %>% 
  

  purrr::set_names('LU01_det_occ_plots',
                   'LU13_det_occ_plots',
                   'LU15_det_occ_plots',
                   'LU21_det_occ_plots') %>% 
  
  purrr::imap(
    ~ggsave(.x,
            file = paste0("figures/",
                          .y,
                          '.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
            dpi = 600,
            width = 12,
            height = 15,
            units = 'in'))

Analysis prep

Read in data

We need the proportional binomial data and the covariate data (from the ACME_camera_script_9-2-2024.R or .Rmd), let’s read those in now and check the structure of each

# response metric (proportional detections from the from the ACME_camera_script_9-2-2024.R or .Rmd)
prop_detections <- read_csv('data/processed/OSM_2022_proportional_detections.csv')
## Rows: 152 Columns: 23
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): site
## dbl (22): black_bear, coyote, fisher, moose, white-tailed_deer, cougar, grey...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# check variable structure
str(prop_detections)
## spc_tbl_ [152 × 23] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ site                    : chr [1:152] "LU01_06" "LU01_10" "LU01_11" "LU01_13" ...
##  $ black_bear              : num [1:152] 7 3 4 7 8 9 4 5 7 7 ...
##  $ coyote                  : num [1:152] 4 4 8 10 11 9 11 0 9 4 ...
##  $ fisher                  : num [1:152] 5 3 3 3 2 1 1 1 0 3 ...
##  $ moose                   : num [1:152] 3 2 5 9 1 0 2 4 1 0 ...
##  $ white-tailed_deer       : num [1:152] 12 5 12 12 13 14 15 9 12 10 ...
##  $ cougar                  : num [1:152] 0 0 1 0 1 0 0 0 0 0 ...
##  $ grey_wolf               : num [1:152] 0 0 2 0 0 0 1 0 0 0 ...
##  $ lynx                    : num [1:152] 0 0 1 0 1 1 0 0 0 2 ...
##  $ red_fox                 : num [1:152] 0 0 2 0 0 0 0 0 4 0 ...
##  $ wolverine               : num [1:152] 0 0 0 0 0 0 0 0 0 0 ...
##  $ caribou                 : num [1:152] 0 0 0 0 0 0 0 0 0 0 ...
##  $ absent_black_bear       : num [1:152] 5 3 8 5 4 3 8 7 5 5 ...
##  $ absent_coyote           : num [1:152] 10 1 6 5 3 5 4 15 6 11 ...
##  $ absent_fisher           : num [1:152] 9 2 11 12 12 13 14 14 15 12 ...
##  $ absent_moose            : num [1:152] 11 3 9 6 13 14 13 11 14 15 ...
##  $ absent_white-tailed_deer: num [1:152] 2 0 2 3 1 0 0 6 3 5 ...
##  $ absent_cougar           : num [1:152] 14 5 13 15 13 14 15 15 15 15 ...
##  $ absent_grey_wolf        : num [1:152] 14 5 12 15 14 14 14 15 15 15 ...
##  $ absent_lynx             : num [1:152] 14 5 13 15 13 13 15 15 15 13 ...
##  $ absent_red_fox          : num [1:152] 14 5 12 15 14 14 15 15 11 15 ...
##  $ absent_wolverine        : num [1:152] 14 5 14 15 14 14 15 15 15 15 ...
##  $ absent_caribou          : num [1:152] 14 5 14 15 14 14 15 15 15 15 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   site = col_character(),
##   ..   black_bear = col_double(),
##   ..   coyote = col_double(),
##   ..   fisher = col_double(),
##   ..   moose = col_double(),
##   ..   `white-tailed_deer` = col_double(),
##   ..   cougar = col_double(),
##   ..   grey_wolf = col_double(),
##   ..   lynx = col_double(),
##   ..   red_fox = col_double(),
##   ..   wolverine = col_double(),
##   ..   caribou = col_double(),
##   ..   absent_black_bear = col_double(),
##   ..   absent_coyote = col_double(),
##   ..   absent_fisher = col_double(),
##   ..   absent_moose = col_double(),
##   ..   `absent_white-tailed_deer` = col_double(),
##   ..   absent_cougar = col_double(),
##   ..   absent_grey_wolf = col_double(),
##   ..   absent_lynx = col_double(),
##   ..   absent_red_fox = col_double(),
##   ..   absent_wolverine = col_double(),
##   ..   absent_caribou = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>
# model covariates (merged HFI and VEG data from the ACME_camera_script_9-2-2024.R or .Rmd)
covariates <- read_csv('data/processed/OSM_2022_covariates.csv',
                       
                       # set the column types to read in correctly
                       col_types = cols(array = col_factor(),
                                        camera = col_factor(),
                                        site = col_factor(),
                                        buff_dist = col_factor(),
                                        .default = col_number()))

# check variable structure
str(covariates)
## spc_tbl_ [3,100 × 77] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ array                       : Factor w/ 4 levels "LU13","LU15",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ camera                      : Factor w/ 96 levels "18","15","03",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ site                        : Factor w/ 155 levels "LU13_18","LU13_15",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ buff_dist                   : Factor w/ 20 levels "250","500","750",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ harvest_area                : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ crop                        : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_aband                  : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_oil                    : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ trail                       : num [1:3100] 0 0 NA 0.5 0 ...
##  $ harvest_area_white_zone     : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ conventional_seismic        : num [1:3100] 0.5 0.5 NA 0.5 1 ...
##  $ pipeline                    : num [1:3100] 0 0.5 NA 0 0 0 0.5 0 0 0 ...
##  $ tame_pasture                : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ rough_pasture               : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ rural_residence             : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ transmission_line           : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_gas                    : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ misc_oil_gas_facility       : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ clearing_unknown            : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ vegetated_edge_roads        : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ road_unimproved             : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ road_gravel_1l              : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ road_gravel_2l              : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ truck_trail                 : num [1:3100] 0.5 0 NA 0 0 0 0 0 0 0 ...
##  $ borrowpits                  : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ sump                        : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ borrowpit_wet               : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ cultivation_abandoned       : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ urban_residence             : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ country_residence           : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ recreation                  : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_other                  : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_bitumen                : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_cased                  : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ road_paved_undiv_2l         : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ road_unclassified           : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ runway                      : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ clearing_wellpad_unconfirmed: num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ facility_unknown            : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ borrowpit_dry               : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ grvl_sand_pit               : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ dugout                      : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ lagoon                      : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ open_pit_mine               : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ low_impact_seismic          : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ surrounding_veg             : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ transfer_station            : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ facility_other              : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ vegetated_edge_railways     : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ fruit_vegetables            : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ residence_clearing          : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ cfo                         : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ landfill                    : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_cleared_not_confirmed  : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ oil_gas_plant               : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ urban_industrial            : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ road_paved_1l               : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ road_paved_undiv_1l         : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ road_winter                 : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_cleared_not_drilled    : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ well_unknown                : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ airp_runway                 : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ reservoir                   : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ campground                  : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ canal                       : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ camp_industrial             : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ rlwy_sgl_track              : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ lc_class20                  : num [1:3100] 0.143 0 0 0 0 ...
##  $ lc_class32                  : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
##  $ lc_class33                  : num [1:3100] 0 0.2 0 0 0 ...
##  $ lc_class34                  : num [1:3100] 0 0.2 0 0 0 ...
##  $ lc_class50                  : num [1:3100] 0.286 0 0.333 0 0 ...
##  $ lc_class110                 : num [1:3100] 0.143 0.2 0 0 0.333 ...
##  $ lc_class120                 : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
##  $ lc_class210                 : num [1:3100] 0.429 0.2 0.333 1 0.667 ...
##  $ lc_class220                 : num [1:3100] 0 0 0 0 0 0.25 0 0 0 0 ...
##  $ lc_class230                 : num [1:3100] 0 0.2 0.333 0 0 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   .default = col_number(),
##   ..   array = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   ..   camera = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   ..   site = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   ..   buff_dist = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   ..   harvest_area = col_number(),
##   ..   crop = col_number(),
##   ..   well_aband = col_number(),
##   ..   well_oil = col_number(),
##   ..   trail = col_number(),
##   ..   harvest_area_white_zone = col_number(),
##   ..   conventional_seismic = col_number(),
##   ..   pipeline = col_number(),
##   ..   tame_pasture = col_number(),
##   ..   rough_pasture = col_number(),
##   ..   rural_residence = col_number(),
##   ..   transmission_line = col_number(),
##   ..   well_gas = col_number(),
##   ..   misc_oil_gas_facility = col_number(),
##   ..   clearing_unknown = col_number(),
##   ..   vegetated_edge_roads = col_number(),
##   ..   road_unimproved = col_number(),
##   ..   road_gravel_1l = col_number(),
##   ..   road_gravel_2l = col_number(),
##   ..   truck_trail = col_number(),
##   ..   borrowpits = col_number(),
##   ..   sump = col_number(),
##   ..   borrowpit_wet = col_number(),
##   ..   cultivation_abandoned = col_number(),
##   ..   urban_residence = col_number(),
##   ..   country_residence = col_number(),
##   ..   recreation = col_number(),
##   ..   well_other = col_number(),
##   ..   well_bitumen = col_number(),
##   ..   well_cased = col_number(),
##   ..   road_paved_undiv_2l = col_number(),
##   ..   road_unclassified = col_number(),
##   ..   runway = col_number(),
##   ..   clearing_wellpad_unconfirmed = col_number(),
##   ..   facility_unknown = col_number(),
##   ..   borrowpit_dry = col_number(),
##   ..   grvl_sand_pit = col_number(),
##   ..   dugout = col_number(),
##   ..   lagoon = col_number(),
##   ..   open_pit_mine = col_number(),
##   ..   low_impact_seismic = col_number(),
##   ..   surrounding_veg = col_number(),
##   ..   transfer_station = col_number(),
##   ..   facility_other = col_number(),
##   ..   vegetated_edge_railways = col_number(),
##   ..   fruit_vegetables = col_number(),
##   ..   residence_clearing = col_number(),
##   ..   cfo = col_number(),
##   ..   landfill = col_number(),
##   ..   well_cleared_not_confirmed = col_number(),
##   ..   oil_gas_plant = col_number(),
##   ..   urban_industrial = col_number(),
##   ..   road_paved_1l = col_number(),
##   ..   road_paved_undiv_1l = col_number(),
##   ..   road_winter = col_number(),
##   ..   well_cleared_not_drilled = col_number(),
##   ..   well_unknown = col_number(),
##   ..   airp_runway = col_number(),
##   ..   reservoir = col_number(),
##   ..   campground = col_number(),
##   ..   canal = col_number(),
##   ..   camp_industrial = col_number(),
##   ..   rlwy_sgl_track = col_number(),
##   ..   lc_class20 = col_number(),
##   ..   lc_class32 = col_number(),
##   ..   lc_class33 = col_number(),
##   ..   lc_class34 = col_number(),
##   ..   lc_class50 = col_number(),
##   ..   lc_class110 = col_number(),
##   ..   lc_class120 = col_number(),
##   ..   lc_class210 = col_number(),
##   ..   lc_class220 = col_number(),
##   ..   lc_class230 = col_number()
##   .. )
##  - attr(*, "problems")=<externalptr>

Format covariates

There are too many covariates to include in the models individually and many of them describe similar HFI features. We can use the info from the README file in this repository which includes detailed descriptions from the ABMI human footprints wall to wall data download website for Year 2021 OR in the relevant_literature folder of this repository (HFI_2021_v1_0_Metadata_Final.pdf).

the current version of this code for the purposes of the 2022-2023 report used a merged dataset from 2021-2022 and 2022-2023 data, howver each year of data the variables were extracted slightly differenty from GIS so final version of this code will include a different formatting process which will likely occur in the ACME_camera_script_9-2-2024.R or .Rmd

covariates_grouped <- covariates %>% 
  
  mutate(borrowpits = rowSums(across(contains('borrowpit'))),
         industrial_sites = camp_industrial + oil_gas_plant + open_pit_mine + 
           rowSums(across(contains('facility'))),
         seismic_lines = rowSums(across(contains('seismic'))),
         wellsites = rowSums(across(contains('well'))),
         roads =  rowSums(across(contains('road'))),
         havest_areas = rowSums(across(contains('harvest'))),
         trails = rowSums(across(contains('trail'))),
         residences = rowSums(across(contains('residence'))),
         pasture = rowSums(across(contains('pasture'))),
         other_transportation_features = runway + airp_runway + rlwy_sgl_track + vegetated_edge_railways,
         crops = crop + fruit_vegetables + cultivation_abandoned,
         water = lagoon + reservoir + dugout + canal,
         .keep = 'unused') %>% 
  
  # remove features we don't need
  select(!c(recreation,
            clearing_unknown,
            cfo,
            grvl_sand_pit,
            transfer_station,
            campground,
            surrounding_veg,
            urban_industrial,
            landfill,
            sump,
            water,
            crops,
            other_transportation_features,
            pasture,
            residences
            )) %>% 
  
  # reorder variables
  relocate(c(pipeline,
             transmission_line,
             borrowpits),
           .after = lc_class230)

# see what's left
names(covariates_grouped)
##  [1] "array"             "camera"            "site"             
##  [4] "buff_dist"         "lc_class20"        "lc_class32"       
##  [7] "lc_class33"        "lc_class34"        "lc_class50"       
## [10] "lc_class110"       "lc_class120"       "lc_class210"      
## [13] "lc_class220"       "lc_class230"       "pipeline"         
## [16] "transmission_line" "borrowpits"        "industrial_sites" 
## [19] "seismic_lines"     "wellsites"         "roads"            
## [22] "havest_areas"      "trails"
# check the structure of new data
str(covariates_grouped)
## tibble [3,100 × 23] (S3: tbl_df/tbl/data.frame)
##  $ array            : Factor w/ 4 levels "LU13","LU15",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ camera           : Factor w/ 96 levels "18","15","03",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ site             : Factor w/ 155 levels "LU13_18","LU13_15",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ buff_dist        : Factor w/ 20 levels "250","500","750",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ lc_class20       : num [1:3100] 0.143 0 0 0 0 ...
##  $ lc_class32       : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
##  $ lc_class33       : num [1:3100] 0 0.2 0 0 0 ...
##  $ lc_class34       : num [1:3100] 0 0.2 0 0 0 ...
##  $ lc_class50       : num [1:3100] 0.286 0 0.333 0 0 ...
##  $ lc_class110      : num [1:3100] 0.143 0.2 0 0 0.333 ...
##  $ lc_class120      : num [1:3100] 0 0 0 0 0 0 0 0 0 0 ...
##  $ lc_class210      : num [1:3100] 0.429 0.2 0.333 1 0.667 ...
##  $ lc_class220      : num [1:3100] 0 0 0 0 0 0.25 0 0 0 0 ...
##  $ lc_class230      : num [1:3100] 0 0.2 0.333 0 0 ...
##  $ pipeline         : num [1:3100] 0 0.5 NA 0 0 0 0.5 0 0 0 ...
##  $ transmission_line: num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ borrowpits       : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ industrial_sites : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ seismic_lines    : num [1:3100] 0.5 0.5 NA 0.5 1 ...
##  $ wellsites        : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ roads            : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ havest_areas     : num [1:3100] 0 0 NA 0 0 0 0 0 0 0 ...
##  $ trails           : num [1:3100] 0.5 0 NA 0.5 0 ...
# check summary of new data
summary(covariates_grouped)
##   array         camera          site        buff_dist      lc_class20     
##  LU13:820   27     :  80   LU13_18:  20   250    : 155   Min.   :0.00000  
##  LU15:780   32     :  80   LU13_15:  20   500    : 155   1st Qu.:0.00000  
##  LU21:720   41     :  80   LU13_03:  20   750    : 155   Median :0.00000  
##  LU01:780   36     :  80   LU13_34:  20   1000   : 155   Mean   :0.02914  
##             16     :  60   LU13_57:  20   1250   : 155   3rd Qu.:0.04545  
##             21     :  60   LU13_16:  20   1500   : 155   Max.   :0.37500  
##             (Other):2660   (Other):2980   (Other):2170                    
##    lc_class32          lc_class33         lc_class34         lc_class50    
##  Min.   :0.000e+00   Min.   :0.000000   Min.   :0.000000   Min.   :0.0000  
##  1st Qu.:0.000e+00   1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.:0.1535  
##  Median :0.000e+00   Median :0.000000   Median :0.004525   Median :0.2500  
##  Mean   :4.373e-05   Mean   :0.012162   Mean   :0.024762   Mean   :0.2733  
##  3rd Qu.:0.000e+00   3rd Qu.:0.009346   3rd Qu.:0.021739   3rd Qu.:0.3750  
##  Max.   :2.273e-02   Max.   :0.400000   Max.   :0.500000   Max.   :0.9091  
##                                                                            
##   lc_class110       lc_class120        lc_class210      lc_class220     
##  Min.   :0.00000   Min.   :0.000000   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:0.04918   1st Qu.:0.000000   1st Qu.:0.1270   1st Qu.:0.05882  
##  Median :0.13333   Median :0.000000   Median :0.1901   Median :0.17903  
##  Mean   :0.13762   Mean   :0.001355   Mean   :0.2120   Mean   :0.16350  
##  3rd Qu.:0.20000   3rd Qu.:0.000000   3rd Qu.:0.2568   3rd Qu.:0.25000  
##  Max.   :0.66667   Max.   :0.200000   Max.   :1.0000   Max.   :0.66667  
##                                                                         
##   lc_class230         pipeline       transmission_line    borrowpits     
##  Min.   :0.00000   Min.   :0.00000   Min.   :0.000000   Min.   :0.00000  
##  1st Qu.:0.06667   1st Qu.:0.00000   1st Qu.:0.000000   1st Qu.:0.00000  
##  Median :0.13208   Median :0.03637   Median :0.000000   Median :0.00000  
##  Mean   :0.14614   Mean   :0.06940   Mean   :0.004712   Mean   :0.01108  
##  3rd Qu.:0.21053   3rd Qu.:0.10638   3rd Qu.:0.000000   3rd Qu.:0.01807  
##  Max.   :0.66667   Max.   :1.00000   Max.   :0.500000   Max.   :0.16667  
##                    NA's   :8         NA's   :8          NA's   :8        
##  industrial_sites   seismic_lines      wellsites           roads        
##  Min.   :0.000000   Min.   :0.0000   Min.   :0.00000   Min.   :0.00000  
##  1st Qu.:0.000000   1st Qu.:0.2774   1st Qu.:0.01541   1st Qu.:0.00000  
##  Median :0.000000   Median :0.3868   Median :0.04408   Median :0.05939  
##  Mean   :0.001448   Mean   :0.4173   Mean   :0.05748   Mean   :0.15189  
##  3rd Qu.:0.000000   3rd Qu.:0.5000   3rd Qu.:0.08125   3rd Qu.:0.27978  
##  Max.   :0.111111   Max.   :1.0000   Max.   :0.50000   Max.   :0.83333  
##  NA's   :8          NA's   :8        NA's   :8         NA's   :8        
##   havest_areas         trails      
##  Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:0.00000   1st Qu.:0.0617  
##  Median :0.00000   Median :0.1522  
##  Mean   :0.04801   Mean   :0.1874  
##  3rd Qu.:0.04506   3rd Qu.:0.2712  
##  Max.   :0.83333   Max.   :1.0000  
##  NA's   :8         NA's   :8
# there are some NAs in the data which will cause problems with modeling/visualization of data ignore for now but will explore these sites specifically after report

covariates_grouped <- covariates_grouped %>% 
  
  # remove rows with NAs
  na.omit()

Subset data & correlation plots

Marissa try to get the purrr code for this to work later

Now we need to subset the data for each buffer width, and then in the same loop let’s make correlation plots for these variables within each buffer

# Couldn't get this to work in purrr yet so using a loop to subset the data, create the plots, and save them all in one section... NEAT

buffer_frames<-list()

for (i in unique(covariates_grouped$buff_dist)){
  
  print(i)
  
  #Subset data based on radius
  df<-covariates_grouped%>%
    filter(buff_dist == i)
  
  #rename dataframe on the fly
  assign(paste("df", i, sep ="_"), df)
  
  #list of dataframes
  buffer_frames<-c(buffer_frames, list(df))
  
  #Subset data based on radius
  df<-covariates_grouped%>%
    filter(buff_dist == i)%>%
    select(where(is.numeric))
  
 #compute a correlation matrix (watch for errors)
 matrix<-cor(df)
 
 #print and save the correlation plot on the go
 #renaming for each buffer as we do
 png(file.path("figures/", paste("correlation_", i, ".png")))
 corrplot::corrplot(matrix,
                    type = 'upper',
                    tl.col = 'black',
                    title = paste0('Variable correlation plot at ', i))
 dev.off()
  
}
## [1] "250"
## Warning in cor(df): the standard deviation is zero
## [1] "500"
## Warning in cor(df): the standard deviation is zero
## [1] "750"
## Warning in cor(df): the standard deviation is zero
## [1] "1000"
## Warning in cor(df): the standard deviation is zero
## [1] "1250"
## Warning in cor(df): the standard deviation is zero
## [1] "1500"
## Warning in cor(df): the standard deviation is zero
## [1] "1750"
## Warning in cor(df): the standard deviation is zero
## [1] "2000"
## Warning in cor(df): the standard deviation is zero
## [1] "2250"
## Warning in cor(df): the standard deviation is zero
## [1] "2500"
## Warning in cor(df): the standard deviation is zero
## [1] "2750"
## Warning in cor(df): the standard deviation is zero
## [1] "3000"
## Warning in cor(df): the standard deviation is zero
## [1] "3250"
## Warning in cor(df): the standard deviation is zero
## [1] "3500"
## Warning in cor(df): the standard deviation is zero
## [1] "3750"
## [1] "4000"
## [1] "4250"
## [1] "4500"
## [1] "4750"
## [1] "5000"
# name list objects so we can extract names for plotting 

buffer_frames <- buffer_frames %>% 
  
  # absurdly long way to do this but for sake of time fuck it
  purrr::set_names('250 meter buffer',
                   '500 meter buffer',
                   '750 meter buffer',
                   '1000 meter buffer',
                   '1250 meter buffer',
                   '1500 meter buffer',
                   '1750 meter buffer',
                   '2000 meter buffer',
                   '2250 meter buffer',
                   '2500 meter buffer',
                   '2750 meter buffer',
                   '3000 meter buffer',
                   '3250 meter buffer',
                   '3500 meter buffer',
                   '3750 meter buffer',
                   '4000 meter buffer',
                   '4250 meter buffer',
                   '4500 meter buffer',
                   '4750 meter buffer',
                   '5000 meter buffer')

Exploratory plots

add more to this section in later when we have more time to explore the covariates and choose which should be inlcuded etc.

hfi_histograms <- buffer_frames %>% 
  
  purrr::imap(
    ~.x %>% 
      
      # filter to just the HFI variables 
      select(where(is.numeric) &
          ! starts_with('lc_class')) %>% 
      
      # pipe into hist.data.frame function to make histograms for each variable
      hist.data.frame(mtitl = paste0('Histograms of HFI variables at ', .y)))

Now let’s do the same thing with the landcover variables

lc_histograms <- buffer_frames %>% 
  
  purrr::imap(
    ~.x %>% 
      
      # filter to just the landcover variables 
      select(where(is.numeric) &
          starts_with('lc_class')) %>% 
      
      # pipe into hist.data.frame function to make histograms for each variable
      hist.data.frame(mtitl = paste0('Histograms of landcover variables at ', .y)))

Add response metric

Now that we have the covariate data formatted we need to add the response metric (monthly proportional presence/absence) to the data frames

final_df <- buffer_frames %>% 
  
  purrr::map(
    ~.x %>% 
      
      left_join(prop_detections,
                by = 'site'))

Analysis

Now we are going to run a global model which includes all HFI and LC variables that at first glance (will do a more thorough check later) seem to have enough data to include as covariates for each buffer width, and then we will compare these models see which buffer width best fit the data for each species.

We don’t need to do ALL the species since many don’t have enough data.

Refer to the ACME_camera_script_9-2-2024.html or .Rmd the plot for proportional monthly detections should provide info on which species we have enough data for, can be found under Response metrics/3.Proportion monthly detections

A brief look at this fig indicates that we have enough for all the mammals in the prop_detections data frame except

  • cougar
  • wolverine
  • caribou??? (may have enough may not)

Black bear

ANDREW/MARISSA FIX later

there is probably a way to shorten the following code to select particular species, I saw Andrew’s for loop in the draft script he wrote but couldn’t quite figure it out so I did this instead, maybe we can merge approaches?

Global model

Let’s start with bears and use purrr to create a global model for every buffer distance

black_bear_mods <- final_df %>%

  # use purrr map to fun the same functions for all buffer sizes ((all objects in list))
  purrr::map(
    ~.x %>%

      # glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
     glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))

Model selection

We will use the model.sel() function from the MuMIn package to compare the global models for each buffer width and see which buffer fits the bear data best

model.sel(black_bear_mods)
## Warning in model.sel.default(black_bear_mods): models are not all fitted to the
## same data
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 250 meter buffer    -0.15620          +   0.6632      0.1236   -2.907000
## 750 meter buffer    -1.95600          +   3.2370     -0.3414    1.344000
## 4250 meter buffer    3.99200          +  -8.6930     -6.5210   -3.613000
## 4500 meter buffer    3.02400          +  -7.7200     -5.4380   -3.311000
## 500 meter buffer    -0.37050          +  -2.7320     -0.6889   -0.670300
## 4000 meter buffer    3.41900          +  -5.5320     -5.3720   -3.507000
## 1000 meter buffer   -0.72430          +   3.9400     -1.4870   -0.003855
## 4750 meter buffer    3.99500          +  -4.9810     -5.4840   -3.628000
## 3500 meter buffer    4.61900          +  -2.5720     -6.4170   -5.508000
## 1250 meter buffer   -0.05382          +  -2.4740     -1.9990   -0.202500
## 5000 meter buffer    4.33400          +  -3.2840     -5.4490   -5.221000
## 3750 meter buffer    4.15700          +  -1.9330     -5.7150   -5.938000
## 2500 meter buffer    3.60400          +   4.0000     -4.3810   -2.807000
## 2000 meter buffer    3.34500          +   2.8120     -5.0190   -2.022000
## 2750 meter buffer    3.32100          +   0.6033     -4.7360   -3.172000
## 1500 meter buffer    2.42300          +   0.1242     -4.4740   -2.707000
## 2250 meter buffer    2.80700          +   4.1350     -4.2460   -2.002000
## 3250 meter buffer    3.42500          +  -2.4720     -5.1480   -3.925000
## 3000 meter buffer    2.95400          +  -2.5420     -4.9970   -3.353000
## 1750 meter buffer    3.07800          +   2.5640     -4.6540   -1.829000
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 250 meter buffer      -0.8100      1.2680     -0.9658     -0.5266    -0.01751
## 750 meter buffer       1.0660      2.4280      1.5200      0.7237     0.72120
## 4250 meter buffer     -4.2940     -5.5240     -5.6240    -32.9500    -5.98100
## 4500 meter buffer     -2.8150     -4.7320     -4.7080    -37.3700    -4.95100
## 500 meter buffer      -0.2043      1.2970     -0.1915     -0.6064    -0.49080
## 4000 meter buffer     -3.4470     -4.8610     -4.7940    -30.8300    -5.16900
## 1000 meter buffer     -0.5946      0.5603     -0.7501      0.1574    -0.69290
## 4750 meter buffer     -3.7660     -4.7590     -5.2050    -44.9600    -5.66400
## 3500 meter buffer     -5.3440     -5.5890     -5.4340    -22.3400    -5.51500
## 1250 meter buffer     -1.6130     -0.1985     -1.7500     -0.6162    -1.33100
## 5000 meter buffer     -4.3580     -4.6980     -5.6180    -50.0300    -5.91900
## 3750 meter buffer     -4.8560     -5.3170     -4.6560    -24.3400    -5.12500
## 2500 meter buffer     -3.4410     -4.7590     -4.3690    -19.7200    -4.84600
## 2000 meter buffer     -3.6160     -4.4010     -4.1740    -11.3800    -4.56800
## 2750 meter buffer     -3.5310     -4.1790     -3.4030    -18.8100    -3.99400
## 1500 meter buffer     -3.6120     -3.1860     -3.7450     -5.4970    -3.89000
## 2250 meter buffer     -2.6530     -4.1160     -3.7100    -14.4200    -4.07700
## 3250 meter buffer     -4.0780     -4.2550     -4.2520    -17.3300    -4.34500
## 3000 meter buffer     -3.6690     -3.7370     -3.8960    -13.8300    -3.79600
## 1750 meter buffer     -4.1040     -3.6840     -4.0900     -8.0540    -4.40300
##                   cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df   logLik
## 250 meter buffer  -0.31540 -0.15940      -0.2193  0.12830  -1.1400 15 -281.417
## 750 meter buffer   0.25910  0.22470       0.2466  0.54550  -0.4503 15 -313.123
## 4250 meter buffer  0.19010  1.32400       1.6020  1.50300  -0.3255 15 -314.209
## 4500 meter buffer  0.54310  1.22800       1.6070  1.48000  -0.9885 15 -315.213
## 500 meter buffer  -0.40100 -0.13060      -0.1239 -0.10060   0.4296 15 -315.684
## 4000 meter buffer -0.06547  0.87580       1.2940  1.32900  -0.6553 15 -315.865
## 1000 meter buffer  0.37460  0.75610       0.7682  0.81240  -1.2510 15 -315.937
## 4750 meter buffer -0.50120  0.45590       0.9784  0.71950   0.1911 15 -316.072
## 3500 meter buffer -1.33600  0.18180       0.7121  0.40060   2.2800 15 -316.302
## 1250 meter buffer  0.90570  1.08600       0.8048  1.25700  -1.4600 15 -316.465
## 5000 meter buffer -0.35550  0.28690       0.8817  0.80080  -0.5400 15 -316.505
## 3750 meter buffer -0.88800  0.14020       0.6929  0.66420   1.2280 15 -316.593
## 2500 meter buffer -0.70050 -0.10240       0.3952  0.34430   0.9719 15 -316.754
## 2000 meter buffer -0.31350  0.13820       0.5285  0.32790   0.5836 15 -316.992
## 2750 meter buffer -0.55870 -0.39300       0.1135  0.07635   1.9580 15 -317.110
## 1500 meter buffer  0.28900  0.71010       0.8328  1.12900  -0.5994 15 -317.112
## 2250 meter buffer -0.40770  0.09159       0.5077  0.49230   0.7172 15 -317.132
## 3250 meter buffer -1.44300  0.10820       0.6210  0.47640   2.1430 15 -317.557
## 3000 meter buffer -1.10100  0.19610       0.6168  0.44360   2.1800 15 -317.912
## 1750 meter buffer -0.33480  0.27400       0.6091  0.44150   0.5543 15 -318.429
##                    AICc delta weight
## 250 meter buffer  596.6  0.00      1
## 750 meter buffer  659.8 63.19      0
## 4250 meter buffer 661.9 65.36      0
## 4500 meter buffer 664.0 67.37      0
## 500 meter buffer  664.9 68.31      0
## 4000 meter buffer 665.3 68.67      0
## 1000 meter buffer 665.4 68.82      0
## 4750 meter buffer 665.7 69.09      0
## 3500 meter buffer 666.1 69.55      0
## 1250 meter buffer 666.5 69.88      0
## 5000 meter buffer 666.5 69.95      0
## 3750 meter buffer 666.7 70.13      0
## 2500 meter buffer 667.0 70.45      0
## 2000 meter buffer 667.5 70.93      0
## 2750 meter buffer 667.7 71.16      0
## 1500 meter buffer 667.8 71.17      0
## 2250 meter buffer 667.8 71.21      0
## 3250 meter buffer 668.6 72.06      0
## 3000 meter buffer 669.4 72.77      0
## 1750 meter buffer 670.4 73.80      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

hmmmm seems fishy to me that the 250 meter buffer which is the only one that had missing data would perform THAT much better than all the others, and really you shouldn’t compare models if they aren’t run on the same data, hence the warning message

Let’s remove the 250 buffer and see what happens

black_bear_mods_no250 <- black_bear_mods  %>% 
  
  # purrr::discard_at will remove an item from a list
  purrr::discard_at('250 meter buffer')


# run model selection again
model.sel(black_bear_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 750 meter buffer    -1.95600          +   3.2370     -0.3414    1.344000
## 4250 meter buffer    3.99200          +  -8.6930     -6.5210   -3.613000
## 4500 meter buffer    3.02400          +  -7.7200     -5.4380   -3.311000
## 500 meter buffer    -0.37050          +  -2.7320     -0.6889   -0.670300
## 4000 meter buffer    3.41900          +  -5.5320     -5.3720   -3.507000
## 1000 meter buffer   -0.72430          +   3.9400     -1.4870   -0.003855
## 4750 meter buffer    3.99500          +  -4.9810     -5.4840   -3.628000
## 3500 meter buffer    4.61900          +  -2.5720     -6.4170   -5.508000
## 1250 meter buffer   -0.05382          +  -2.4740     -1.9990   -0.202500
## 5000 meter buffer    4.33400          +  -3.2840     -5.4490   -5.221000
## 3750 meter buffer    4.15700          +  -1.9330     -5.7150   -5.938000
## 2500 meter buffer    3.60400          +   4.0000     -4.3810   -2.807000
## 2000 meter buffer    3.34500          +   2.8120     -5.0190   -2.022000
## 2750 meter buffer    3.32100          +   0.6033     -4.7360   -3.172000
## 1500 meter buffer    2.42300          +   0.1242     -4.4740   -2.707000
## 2250 meter buffer    2.80700          +   4.1350     -4.2460   -2.002000
## 3250 meter buffer    3.42500          +  -2.4720     -5.1480   -3.925000
## 3000 meter buffer    2.95400          +  -2.5420     -4.9970   -3.353000
## 1750 meter buffer    3.07800          +   2.5640     -4.6540   -1.829000
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 750 meter buffer       1.0660      2.4280      1.5200      0.7237      0.7212
## 4250 meter buffer     -4.2940     -5.5240     -5.6240    -32.9500     -5.9810
## 4500 meter buffer     -2.8150     -4.7320     -4.7080    -37.3700     -4.9510
## 500 meter buffer      -0.2043      1.2970     -0.1915     -0.6064     -0.4908
## 4000 meter buffer     -3.4470     -4.8610     -4.7940    -30.8300     -5.1690
## 1000 meter buffer     -0.5946      0.5603     -0.7501      0.1574     -0.6929
## 4750 meter buffer     -3.7660     -4.7590     -5.2050    -44.9600     -5.6640
## 3500 meter buffer     -5.3440     -5.5890     -5.4340    -22.3400     -5.5150
## 1250 meter buffer     -1.6130     -0.1985     -1.7500     -0.6162     -1.3310
## 5000 meter buffer     -4.3580     -4.6980     -5.6180    -50.0300     -5.9190
## 3750 meter buffer     -4.8560     -5.3170     -4.6560    -24.3400     -5.1250
## 2500 meter buffer     -3.4410     -4.7590     -4.3690    -19.7200     -4.8460
## 2000 meter buffer     -3.6160     -4.4010     -4.1740    -11.3800     -4.5680
## 2750 meter buffer     -3.5310     -4.1790     -3.4030    -18.8100     -3.9940
## 1500 meter buffer     -3.6120     -3.1860     -3.7450     -5.4970     -3.8900
## 2250 meter buffer     -2.6530     -4.1160     -3.7100    -14.4200     -4.0770
## 3250 meter buffer     -4.0780     -4.2550     -4.2520    -17.3300     -4.3450
## 3000 meter buffer     -3.6690     -3.7370     -3.8960    -13.8300     -3.7960
## 1750 meter buffer     -4.1040     -3.6840     -4.0900     -8.0540     -4.4030
##                   cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df   logLik
## 750 meter buffer   0.25910  0.22470       0.2466  0.54550  -0.4503 15 -313.123
## 4250 meter buffer  0.19010  1.32400       1.6020  1.50300  -0.3255 15 -314.209
## 4500 meter buffer  0.54310  1.22800       1.6070  1.48000  -0.9885 15 -315.213
## 500 meter buffer  -0.40100 -0.13060      -0.1239 -0.10060   0.4296 15 -315.684
## 4000 meter buffer -0.06547  0.87580       1.2940  1.32900  -0.6553 15 -315.865
## 1000 meter buffer  0.37460  0.75610       0.7682  0.81240  -1.2510 15 -315.937
## 4750 meter buffer -0.50120  0.45590       0.9784  0.71950   0.1911 15 -316.072
## 3500 meter buffer -1.33600  0.18180       0.7121  0.40060   2.2800 15 -316.302
## 1250 meter buffer  0.90570  1.08600       0.8048  1.25700  -1.4600 15 -316.465
## 5000 meter buffer -0.35550  0.28690       0.8817  0.80080  -0.5400 15 -316.505
## 3750 meter buffer -0.88800  0.14020       0.6929  0.66420   1.2280 15 -316.593
## 2500 meter buffer -0.70050 -0.10240       0.3952  0.34430   0.9719 15 -316.754
## 2000 meter buffer -0.31350  0.13820       0.5285  0.32790   0.5836 15 -316.992
## 2750 meter buffer -0.55870 -0.39300       0.1135  0.07635   1.9580 15 -317.110
## 1500 meter buffer  0.28900  0.71010       0.8328  1.12900  -0.5994 15 -317.112
## 2250 meter buffer -0.40770  0.09159       0.5077  0.49230   0.7172 15 -317.132
## 3250 meter buffer -1.44300  0.10820       0.6210  0.47640   2.1430 15 -317.557
## 3000 meter buffer -1.10100  0.19610       0.6168  0.44360   2.1800 15 -317.912
## 1750 meter buffer -0.33480  0.27400       0.6091  0.44150   0.5543 15 -318.429
##                    AICc delta weight
## 750 meter buffer  659.8  0.00  0.504
## 4250 meter buffer 661.9  2.17  0.170
## 4500 meter buffer 664.0  4.18  0.062
## 500 meter buffer  664.9  5.12  0.039
## 4000 meter buffer 665.3  5.48  0.032
## 1000 meter buffer 665.4  5.63  0.030
## 4750 meter buffer 665.7  5.90  0.026
## 3500 meter buffer 666.1  6.36  0.021
## 1250 meter buffer 666.5  6.68  0.018
## 5000 meter buffer 666.5  6.76  0.017
## 3750 meter buffer 666.7  6.94  0.016
## 2500 meter buffer 667.0  7.26  0.013
## 2000 meter buffer 667.5  7.74  0.011
## 2750 meter buffer 667.7  7.97  0.009
## 1500 meter buffer 667.8  7.98  0.009
## 2250 meter buffer 667.8  8.02  0.009
## 3250 meter buffer 668.6  8.87  0.006
## 3000 meter buffer 669.4  9.58  0.004
## 1750 meter buffer 670.4 10.61  0.003
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

this looks much more realistic; the 500 m buffer is top model for black bears

So what we will do for each species is remove the 250 meter buffer for now since there are some data missing. and compare just the other buffer sizes that contain the full data set

Model summary

Let’s take a look at the model summary for the top model

summary(black_bear_mods_no250$`500 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(black_bear, absent_black_bear) ~ seismic_lines + pipeline +  
##     borrowpits + wellsites + roads + trails + lc_class20 + lc_class34 +  
##     lc_class50 + lc_class110 + lc_class210 + lc_class220 + lc_class230 +  
##     (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    661.4    706.7   -315.7    631.4      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.03008  0.1735  
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -0.3705     1.6088  -0.230    0.818
## seismic_lines  -0.1239     0.3770  -0.329    0.742
## pipeline       -0.4010     0.5776  -0.694    0.488
## borrowpits     -2.7317     3.5769  -0.764    0.445
## wellsites       0.4296     0.7900   0.544    0.587
## roads          -0.1306     0.5173  -0.252    0.801
## trails         -0.1006     0.4442  -0.226    0.821
## lc_class20     -0.6703     1.8809  -0.356    0.722
## lc_class34     -0.6064     1.7885  -0.339    0.735
## lc_class50     -0.4908     1.5642  -0.314    0.754
## lc_class110    -0.6889     1.5636  -0.441    0.660
## lc_class210    -0.2043     1.5477  -0.132    0.895
## lc_class220     1.2969     1.5961   0.813    0.416
## lc_class230    -0.1915     1.5672  -0.122    0.903

Let’s repeat this process for each species that we have enough data for.

Caribou

We may or may not have enough data for caribou but let’s try it at least for this preliminary report

We can use the same code from black beasr (above) to run global models for each buffer width except remember we want to remove 250 meters

And in the same chunk to save time let’s also run the model.sel() function

caribou_mods_no250 <- final_df %>%
  
  # remove 350 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(caribou, absent_caribou) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
# model selection
model.sel(caribou_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 4250 meter buffer    -27.510          +   -2.572      10.680      23.580
## 5000 meter buffer    -26.110          +  -14.020      11.220      28.530
## 3250 meter buffer    -54.520          +    5.926      45.440      48.940
## 3750 meter buffer    -22.970          +    1.763      18.410      18.890
## 2500 meter buffer    -55.780          +  -17.250      49.990      43.490
## 4000 meter buffer    -23.790          +   -2.363      12.090      18.610
## 1000 meter buffer    -49.330          +    1.913      38.170      38.320
## 2250 meter buffer    -37.490          +   10.600      29.720      28.270
## 1250 meter buffer    -35.560          +  -17.040      28.720      29.010
## 750 meter buffer     -32.650          +    7.003      21.360      31.990
## 2000 meter buffer    -46.300          +   -7.312      37.500      38.610
## 1500 meter buffer    -35.700          +   -1.047      26.290      32.700
## 1750 meter buffer    -38.510          +  -10.180      33.360      29.810
## 2750 meter buffer    -37.530          +    4.207      30.330      29.170
## 500 meter buffer      -7.754          +   16.820      -2.579       5.773
## 3000 meter buffer    -45.530          +    2.298      37.320      38.900
## 3500 meter buffer    -63.670          +  -13.310      54.360      59.490
## 4500 meter buffer    -21.800          +   13.100       7.141      18.740
## 4750 meter buffer    -25.250          +   -9.999       8.156      28.600
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 4250 meter buffer      13.220     -6.7470      9.5450     78.5300     13.4800
## 5000 meter buffer      16.220     -5.4270     17.5100     73.8200     16.2600
## 3250 meter buffer      46.600     27.9800     47.0200     76.6500     47.4500
## 3750 meter buffer      19.020     -3.2270     11.3100     31.1700     15.4400
## 2500 meter buffer      52.550     37.0000     42.0200     70.4000     52.3200
## 4000 meter buffer      14.750     -4.9520      8.2450     64.7000     12.7800
## 1000 meter buffer      46.540     36.1700     31.9900     33.9900     40.6700
## 2250 meter buffer      34.300     18.0700     25.1000     42.9400     29.5900
## 1250 meter buffer      34.540     26.7800     23.4700     31.5500     31.6300
## 750 meter buffer       27.740     26.4000     18.3100     20.8200     25.6600
## 2000 meter buffer      42.080     29.6500     33.9100     48.4000     39.6300
## 1500 meter buffer      29.300     26.5500     28.3300     39.7500     32.3900
## 1750 meter buffer      34.380     26.0100     27.4400     38.6000     34.6200
## 2750 meter buffer      32.490     20.6700     28.2200     60.1800     32.8600
## 500 meter buffer        1.668     -0.3208     -0.2206      0.5675      0.7308
## 3000 meter buffer      39.700     23.6100     39.3800     70.7300     39.0800
## 3500 meter buffer      55.940     34.6700     51.7100     74.0700     55.0100
## 4500 meter buffer       6.660     -7.1570      5.4240     58.7800      9.9760
## 4750 meter buffer      12.710     -6.9930     11.5200     74.1400     12.8000
##                    cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df  logLik
## 4250 meter buffer   -14.980   19.870     14.09000   13.890   29.300 15 -65.472
## 5000 meter buffer    -7.281   14.060     10.08000    8.032   26.300 15 -69.010
## 3250 meter buffer    -4.064    9.023      5.37200    7.514   10.280 15 -69.014
## 3750 meter buffer    -5.339    8.001      6.09400    8.063    3.877 15 -69.515
## 2500 meter buffer   -15.920    5.627      1.52100    4.422   14.000 15 -69.588
## 4000 meter buffer    -6.043   14.070     10.27000   12.300   11.120 15 -71.099
## 1000 meter buffer   -13.330    7.164      3.25400    4.758   11.760 15 -72.230
## 2250 meter buffer    -9.221    6.348      4.43300    7.423    1.271 15 -72.657
## 1250 meter buffer   -18.310    1.743     -1.00500    2.149    7.739 15 -73.057
## 750 meter buffer   -673.800    4.845      3.45500    2.515    2.598 15 -73.902
## 2000 meter buffer    -6.345    6.081      2.94200    7.057    3.301 15 -74.197
## 1500 meter buffer   -13.130    1.353     -0.80090    3.410    9.249 15 -74.866
## 1750 meter buffer    -6.629    2.111      0.09778    4.664    2.160 15 -74.960
## 2750 meter buffer    -4.771    4.033      2.13800    6.000    1.521 15 -75.182
## 500 meter buffer  -1084.000    2.307      2.26100    4.144   -1.979 15 -88.060
## 3000 meter buffer    -6.259    7.011      4.02500    6.854    5.822 15        
## 3500 meter buffer    -1.937   12.200      7.68700    8.994   13.840 15        
## 4500 meter buffer   -30.210   17.220     12.90000    9.281   39.010 15        
## 4750 meter buffer   -18.350   17.540     12.63000    9.146   41.030 15        
##                    AICc delta
## 4250 meter buffer 164.5  0.00
## 5000 meter buffer 171.5  7.08
## 3250 meter buffer 171.6  7.08
## 3750 meter buffer 172.6  8.08
## 2500 meter buffer 172.7  8.23
## 4000 meter buffer 175.7 11.25
## 1000 meter buffer 178.0 13.52
## 2250 meter buffer 178.8 14.37
## 1250 meter buffer 179.6 15.17
## 750 meter buffer  181.3 16.86
## 2000 meter buffer 181.9 17.45
## 1500 meter buffer 183.3 18.79
## 1750 meter buffer 183.4 18.97
## 2750 meter buffer 183.9 19.42
## 500 meter buffer  209.6 45.17
## 3000 meter buffer            
## 3500 meter buffer            
## 4500 meter buffer            
## 4750 meter buffer            
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

We get a warning that there are some model convergence problems, I expect this is because we don’t have enough data for caribou but I don’t have time to dig into this now so we will investigate more closely for final analysis

For caribou 1250m buffer is top model for now

Let’s take a closer look at the top model summary

summary(caribou_mods_no250$`1250 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(caribou, absent_caribou) ~ seismic_lines + pipeline + borrowpits +  
##     wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +  
##     lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 |      array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    176.1    221.5    -73.1    146.1      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 7.538    2.746   
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)    -35.563     11.737  -3.030  0.00245 **
## seismic_lines   -1.005      1.935  -0.519  0.60351   
## pipeline       -18.310      7.026  -2.606  0.00916 **
## borrowpits     -17.036     13.559  -1.256  0.20894   
## wellsites        7.739      4.257   1.818  0.06911 . 
## roads            1.743      2.439   0.715  0.47485   
## trails           2.149      2.523   0.852  0.39436   
## lc_class20      29.014     12.595   2.304  0.02124 * 
## lc_class34      31.546     11.174   2.823  0.00475 **
## lc_class50      31.628     11.987   2.639  0.00833 **
## lc_class110     28.725     11.951   2.404  0.01624 * 
## lc_class210     34.536     12.811   2.696  0.00702 **
## lc_class220     26.785     11.672   2.295  0.02174 * 
## lc_class230     23.469     11.507   2.040  0.04139 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There’s nothing that catches my eye immediately as being sus about this particular model so it may not have been one with convergence issues. We will keep it in report for now

Coyote

coyote_mods_no250 <- final_df %>%
  
  # remove 350 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(coyote, absent_coyote) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(coyote_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int))  cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 4750 meter buffer   -0.08495          + -10.16000     -0.2466       5.790
## 5000 meter buffer   -0.21630          +  -8.75800      0.1736       5.630
## 4500 meter buffer   -1.47500          + -12.23000      0.0884       7.496
## 4250 meter buffer   -2.21600          +  -7.04100      0.9704       8.918
## 3000 meter buffer   -3.11000          + -10.63000      2.4960       6.791
## 2750 meter buffer   -1.59100          +  -8.08900      1.2310       6.159
## 3500 meter buffer   -2.12900          +  -5.66000      1.3850       6.319
## 3250 meter buffer   -3.29900          +  -8.29300      2.3430       7.866
## 4000 meter buffer   -2.76500          +  -0.01952      2.1350       8.384
## 2500 meter buffer   -3.37600          +  -8.08200      2.5140       6.989
## 3750 meter buffer   -3.07500          +  -2.89100      2.3580       6.734
## 2250 meter buffer   -2.43700          +  -3.13700      1.6200       5.359
## 1750 meter buffer   -3.05600          +  -4.40500      0.7562       6.774
## 1500 meter buffer   -1.16400          +  -5.16200     -0.9758       3.042
## 2000 meter buffer   -1.69600          +  -2.50300     -0.2595       4.830
## 1250 meter buffer   -1.50700          +  -1.89200     -0.3376       3.273
## 500 meter buffer     3.04700          + -11.87000     -4.1530      -3.354
## 1000 meter buffer   -0.92770          +  -2.30100     -0.6898       2.190
## 750 meter buffer     0.80740          +   0.55220     -2.3420      -3.307
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 4750 meter buffer      5.9710     2.89600      -6.603      8.9250     -1.2880
## 5000 meter buffer      6.6510     3.32100      -6.293      7.9440     -0.7587
## 4500 meter buffer      6.8050     3.39300      -5.209     13.2700     -0.2868
## 4250 meter buffer      6.0770     1.84200      -3.960     -5.1260     -0.5428
## 3000 meter buffer      7.2650     4.32100      -5.319     -4.7030      0.8320
## 2750 meter buffer      4.4880     1.25100      -5.400    -12.3500     -0.9980
## 3500 meter buffer      5.3470     0.93400      -5.181    -10.4200     -0.9148
## 3250 meter buffer      6.7610     2.39800      -4.559     -8.6000      0.2598
## 4000 meter buffer      6.0260     0.63670      -2.605    -32.8000     -0.8201
## 2500 meter buffer      5.3330     2.92300      -4.418     -6.2550      0.4718
## 3750 meter buffer      6.1540     1.38400      -3.555    -14.0500     -0.1334
## 2250 meter buffer      4.1120     1.09100      -3.254     -7.7280     -0.3658
## 1750 meter buffer      1.9430     1.53400      -1.629     -0.6329     -0.3105
## 1500 meter buffer      0.3031    -0.95300      -3.445     -2.6330     -2.9080
## 2000 meter buffer      2.0460     0.04994      -3.122     -4.0260     -1.3090
## 1250 meter buffer     -0.1270    -0.35720      -2.318     -2.9060     -1.8130
## 500 meter buffer      -4.0790    -4.09700      -3.968     -3.3830     -4.7730
## 1000 meter buffer     -0.5428    -0.63960      -1.230     -0.5333     -1.6720
## 750 meter buffer      -2.2730    -1.87100      -3.069     -0.4518     -2.6920
##                   cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df   logLik
## 4750 meter buffer  2.05600 -1.81500    -0.999200 -6.62800   4.1390 15 -299.115
## 5000 meter buffer  1.38100 -2.21200    -1.348000 -7.12300   4.5840 15 -300.974
## 4500 meter buffer  2.62000 -0.94760    -0.540000 -5.50800   4.0930 15 -305.186
## 4250 meter buffer -0.03671  0.72540     0.731600 -4.59300   6.3700 15 -309.746
## 3000 meter buffer  2.47400  0.45800    -0.003364 -2.42000   2.6550 15 -312.922
## 2750 meter buffer  1.93500  0.90430     0.451700 -2.22000   3.2820 15 -313.881
## 3500 meter buffer  0.21890  1.39600     1.255000 -2.75000   5.1700 15 -314.162
## 3250 meter buffer  1.40200  1.60100     1.093000 -1.93800   3.0620 15 -314.444
## 4000 meter buffer -0.77580  1.60200     1.544000 -3.12100   5.6220 15 -314.804
## 2500 meter buffer  2.83200  1.73400     0.936700 -1.11700   0.5833 15 -315.340
## 3750 meter buffer -0.08756  1.61400     1.363000 -2.78700   5.5600 15 -315.561
## 2250 meter buffer  2.34200  1.55500     0.536100 -1.01400   0.9714 15 -321.930
## 1750 meter buffer  3.40000  2.40400     0.993400 -0.08169   0.8029 15 -323.037
## 1500 meter buffer  2.45200  2.80400     1.233000  0.15170   1.9900 15 -325.175
## 2000 meter buffer  2.46400  1.82100     0.691800 -0.65260   1.9360 15 -327.232
## 1250 meter buffer  2.25200  2.12100     1.100000  0.06838   1.4950 15 -331.548
## 500 meter buffer   0.41960  0.09539    -0.693400  0.04604   1.3450 15 -336.563
## 1000 meter buffer  1.86600  1.38800     0.365400 -0.41380   0.1284 15 -340.370
## 750 meter buffer   1.09100  0.83070    -0.295400  0.04470   1.6400 15 -341.549
##                    AICc delta weight
## 4750 meter buffer 631.8  0.00  0.863
## 5000 meter buffer 635.5  3.72  0.135
## 4500 meter buffer 643.9 12.14  0.002
## 4250 meter buffer 653.0 21.26  0.000
## 3000 meter buffer 659.4 27.61  0.000
## 2750 meter buffer 661.3 29.53  0.000
## 3500 meter buffer 661.9 30.09  0.000
## 3250 meter buffer 662.4 30.66  0.000
## 4000 meter buffer 663.1 31.38  0.000
## 2500 meter buffer 664.2 32.45  0.000
## 3750 meter buffer 664.7 32.89  0.000
## 2250 meter buffer 677.4 45.63  0.000
## 1750 meter buffer 679.6 47.84  0.000
## 1500 meter buffer 683.9 52.12  0.000
## 2000 meter buffer 688.0 56.23  0.000
## 1250 meter buffer 696.6 64.87  0.000
## 500 meter buffer  706.7 74.90  0.000
## 1000 meter buffer 714.3 82.51  0.000
## 750 meter buffer  716.6 84.87  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

for coyote top model appears to be 4500 m by quite a bit

Let’s get the model summary for this model

summary(coyote_mods_no250$`4500 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(coyote, absent_coyote) ~ seismic_lines + pipeline + borrowpits +  
##     wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +  
##     lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 |      array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    640.4    685.7   -305.2    610.4      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  array  (Intercept) 3.067e-10 1.751e-05
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.4754     2.5310  -0.583   0.5599    
## seismic_lines  -0.5400     1.0561  -0.511   0.6091    
## pipeline        2.6202     1.6554   1.583   0.1135    
## borrowpits    -12.2269     6.0760  -2.012   0.0442 *  
## wellsites       4.0932     2.7281   1.500   0.1335    
## roads          -0.9476     1.1313  -0.838   0.4022    
## trails         -5.5084     1.1641  -4.732 2.22e-06 ***
## lc_class20      7.4961     3.6399   2.059   0.0395 *  
## lc_class34     13.2692    14.4421   0.919   0.3582    
## lc_class50     -0.2868     2.4402  -0.118   0.9064    
## lc_class110     0.0884     2.4431   0.036   0.9711    
## lc_class210     6.8048     2.7885   2.440   0.0147 *  
## lc_class220     3.3933     2.3817   1.425   0.1542    
## lc_class230    -5.2090     2.5779  -2.021   0.0433 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is one lc class with a very high estimate and SE which seems a bit sus to me

Fisher

fisher_mods_no250 <- final_df %>%
  
  # remove 350 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(fisher, absent_fisher) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')
# model selection
model.sel(fisher_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 2500 meter buffer      5.547          +  -22.340    -5.56900      -6.332
## 2750 meter buffer      7.279          +  -18.600    -7.12700      -6.177
## 2000 meter buffer      2.772          +  -20.520    -3.01400      -3.670
## 3250 meter buffer      6.896          +  -19.870    -7.10600      -4.005
## 2250 meter buffer      3.923          +  -20.000    -4.24300      -3.839
## 3000 meter buffer      6.438          +  -25.000    -7.37800      -4.165
## 1750 meter buffer      3.837          +  -17.190    -4.64200      -5.103
## 3750 meter buffer      6.808          +  -23.750    -6.71700      -4.844
## 1500 meter buffer      1.145          +  -14.140    -1.38900      -2.486
## 500 meter buffer    -206.300          +   -3.105   203.20000     204.100
## 3500 meter buffer      7.742          +  -15.450    -7.39800      -3.618
## 4000 meter buffer      4.646          +  -21.930    -4.95100      -1.906
## 750 meter buffer     -11.960          +  -10.280     9.77700      10.670
## 4250 meter buffer      4.041          +  -20.250    -3.99900      -3.637
## 1250 meter buffer     -0.424          +  -10.130    -0.07718      -1.525
## 4750 meter buffer      4.608          +  -12.620    -3.71200      -6.128
## 4500 meter buffer      3.998          +  -14.550    -3.68300      -6.425
## 5000 meter buffer      5.774          +   -9.895    -4.31000      -6.296
## 1000 meter buffer     -2.659          +   -5.299     0.54330       2.036
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 2500 meter buffer     -7.9120     -8.0340      -4.457     -4.5200     -5.8080
## 2750 meter buffer     -9.6060    -10.1400      -4.994    -15.8800     -7.0410
## 2000 meter buffer     -5.4840     -4.1790      -5.517     -3.6230     -2.5970
## 3250 meter buffer    -10.6300     -9.1580      -4.190    -19.9300     -7.8010
## 2250 meter buffer     -7.1950     -5.1290      -4.658     -2.8610     -3.9040
## 3000 meter buffer    -10.3000     -9.1890      -5.101    -11.5900     -7.5710
## 1750 meter buffer     -6.7840     -5.8860      -4.930     -4.9150     -3.9730
## 3750 meter buffer     -8.9080     -9.7550      -5.344    -25.5400     -7.9620
## 1500 meter buffer     -4.3270     -3.5170      -3.169     -3.1480     -1.4970
## 500 meter buffer     202.8000    204.1000     202.100    204.3000    203.1000
## 3500 meter buffer     -9.9350     -9.5030      -5.713    -28.4700     -8.8960
## 4000 meter buffer     -5.5530     -8.6370      -3.874    -26.3800     -6.1570
## 750 meter buffer      10.1900      9.3110       9.756      9.9900      9.4750
## 4250 meter buffer     -4.4600     -8.1490      -3.317    -22.8000     -4.9060
## 1250 meter buffer     -2.7670     -2.5670      -1.689     -0.8512     -0.1994
## 4750 meter buffer     -5.7610     -7.7440      -2.656    -26.0200     -5.4220
## 4500 meter buffer     -4.1480     -8.2130      -3.430    -19.3000     -4.5210
## 5000 meter buffer     -6.4240     -8.5620      -3.544    -42.4300     -6.4730
## 1000 meter buffer     -0.2956     -0.7572      -0.178      1.7690      1.0030
##                   cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl)  cnd(wll) df   logLik
## 2500 meter buffer   0.5878  -3.0710       -3.646  -2.0590   0.93200 15 -168.507
## 2750 meter buffer   0.5857  -3.6490       -4.003  -2.3110   0.41190 15 -168.639
## 2000 meter buffer   2.0200  -2.4280       -2.964  -2.1210   0.55590 15 -168.928
## 3250 meter buffer   1.6390  -3.2260       -3.290  -1.3580  -2.33000 15 -170.390
## 2250 meter buffer   1.4420  -2.7610       -3.150  -2.4170   0.20680 15 -170.610
## 3000 meter buffer   1.3450  -2.1540       -2.697  -0.7910  -1.74200 15 -170.743
## 1750 meter buffer   0.5733  -2.4750       -2.912  -2.2060   3.06800 15 -170.803
## 3750 meter buffer   1.9710  -2.8770       -3.187  -0.9760  -2.16000 15 -171.232
## 1500 meter buffer   1.5190  -1.9160       -2.835  -0.6800  -0.38500 15 -171.327
## 500 meter buffer    1.3750   0.1708        0.206  -0.3670  -0.19210 15 -171.549
## 3500 meter buffer   1.9290  -3.5470       -3.516  -1.9240  -3.04500 15 -172.047
## 4000 meter buffer   3.1890  -2.6640       -2.601  -0.8963  -3.51000 15 -172.099
## 750 meter buffer    0.9276  -1.3280       -1.354  -0.8490   1.64000 15 -172.508
## 4250 meter buffer   3.5880  -3.0850       -3.155  -0.6601  -4.18600 15 -172.554
## 1250 meter buffer   0.9799  -1.7590       -2.175  -0.5384   0.02037 15 -173.142
## 4750 meter buffer   6.8560  -3.8530       -3.420  -0.3058 -11.76000 15 -174.034
## 4500 meter buffer   2.9530  -3.2190       -3.325  -0.7649  -3.86300 15 -174.210
## 5000 meter buffer   3.7570  -4.4220       -3.850  -1.5630  -5.56700 15 -174.883
## 1000 meter buffer   0.4348  -0.9629       -1.086  -0.5941   1.69100 15 -175.884
##                    AICc delta weight
## 2500 meter buffer 370.5  0.00  0.303
## 2750 meter buffer 370.8  0.26  0.266
## 2000 meter buffer 371.4  0.84  0.199
## 3250 meter buffer 374.3  3.77  0.046
## 2250 meter buffer 374.8  4.21  0.037
## 3000 meter buffer 375.0  4.47  0.032
## 1750 meter buffer 375.1  4.59  0.030
## 3750 meter buffer 376.0  5.45  0.020
## 1500 meter buffer 376.2  5.64  0.018
## 500 meter buffer  376.6  6.08  0.014
## 3500 meter buffer 377.6  7.08  0.009
## 4000 meter buffer 377.7  7.18  0.008
## 750 meter buffer  378.5  8.00  0.006
## 4250 meter buffer 378.6  8.09  0.005
## 1250 meter buffer 379.8  9.27  0.003
## 4750 meter buffer 381.6 11.05  0.001
## 4500 meter buffer 381.9 11.41  0.001
## 5000 meter buffer 383.3 12.75  0.001
## 1000 meter buffer 385.3 14.75  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For fisher top model is 2000 meter

Let’s print the summary for this model

summary(fisher_mods_no250$`2000 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(fisher, absent_fisher) ~ seismic_lines + pipeline + borrowpits +  
##     wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +  
##     lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 |      array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    367.9    413.2   -168.9    337.9      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 1.783    1.335   
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)     2.7721     2.5511   1.087   0.2772   
## seismic_lines  -2.9635     1.0988  -2.697   0.0070 **
## pipeline        2.0200     1.4830   1.362   0.1732   
## borrowpits    -20.5174     9.3253  -2.200   0.0278 * 
## wellsites       0.5559     2.5248   0.220   0.8257   
## roads          -2.4279     1.4090  -1.723   0.0849 . 
## trails         -2.1214     1.1269  -1.883   0.0598 . 
## lc_class20     -3.6696     2.8889  -1.270   0.2040   
## lc_class34     -3.6234     7.2996  -0.496   0.6196   
## lc_class50     -2.5967     2.2342  -1.162   0.2451   
## lc_class110    -3.0135     2.4184  -1.246   0.2127   
## lc_class210    -5.4838     2.8300  -1.938   0.0527 . 
## lc_class220    -4.1786     2.4971  -1.673   0.0943 . 
## lc_class230    -5.5168     2.5162  -2.192   0.0283 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Again lc_class34 has a very high standard error, we may not have enough data in this landcover class to use in the final analysis

Grey wolf

wolf_mods_no250 <- final_df %>%
  
  # remove 350 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(grey_wolf, absent_grey_wolf) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(wolf_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 3750 meter buffer    -10.510          +  18.8400      12.540       7.429
## 4000 meter buffer    -10.300          +  19.4600       9.827       7.102
## 3250 meter buffer    -14.660          +  22.6900      16.360       7.025
## 4500 meter buffer     -6.936          +  12.2600       5.408       9.396
## 4250 meter buffer     -5.790          +  15.0200       5.380       6.754
## 3000 meter buffer    -15.640          +  18.8300      16.430       7.908
## 3500 meter buffer    -12.600          +  22.6800      15.560       6.655
## 2750 meter buffer    -13.670          +  24.7600      14.190       7.631
## 4750 meter buffer     -9.296          +   4.2280       6.038       9.575
## 5000 meter buffer     -9.007          +   0.6744       6.215       7.817
## 2500 meter buffer    -11.530          +  17.8500      12.130       7.462
## 1750 meter buffer    -13.220          +  10.1900      11.250       9.449
## 1250 meter buffer    -11.180          +   6.6470       8.402       9.471
## 1500 meter buffer    -10.410          +   6.1110       7.866       8.383
## 2000 meter buffer    -12.200          +  10.8000      10.250      10.510
## 2250 meter buffer    -12.430          +  14.9400      11.980       9.630
## 1000 meter buffer     -7.422          +  19.3300       4.612       6.292
## 750 meter buffer      -5.668          +   8.0980       2.376       2.684
## 500 meter buffer      -5.779          +  -9.0020       3.057       2.157
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 3750 meter buffer      13.000       5.047       6.369     -75.950       8.054
## 4000 meter buffer      13.090       4.363       5.602     -70.960       7.448
## 3250 meter buffer      16.790       9.845      10.600     -47.540      12.100
## 4500 meter buffer       8.556       2.473       3.475     -83.970       4.251
## 4250 meter buffer       7.508       1.306       2.788     -82.100       3.687
## 3000 meter buffer      16.710       9.927      12.690     -35.210      12.860
## 3500 meter buffer      14.910       8.070      10.060     -58.820      10.670
## 2750 meter buffer      14.040       8.286      11.330     -26.540      10.980
## 4750 meter buffer       8.955       4.202       4.422     -68.030       5.646
## 5000 meter buffer       9.197       2.776       4.305     -80.160       5.340
## 2500 meter buffer      10.680       6.564       8.132     -19.930       8.321
## 1750 meter buffer      12.050       6.682       7.176      -6.558       8.335
## 1250 meter buffer       9.470       5.231       5.503      -2.647       5.362
## 1500 meter buffer       8.308       4.247       4.828      -3.690       5.024
## 2000 meter buffer      11.420       6.183       6.063     -11.320       7.652
## 2250 meter buffer      11.170       7.697       8.091     -14.060       8.957
## 1000 meter buffer       4.182       3.625       3.463      -1.872       3.417
## 750 meter buffer        2.134       1.526       1.817      -1.649       1.567
## 500 meter buffer        2.405       2.151       1.530       1.168       1.817
##                    cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df   logLik
## 3750 meter buffer -15.09000  0.10330     -0.09271  -2.2610  15.1000 15 -174.803
## 4000 meter buffer -15.41000  0.97660      0.60560  -1.1840  13.8500 15 -176.036
## 3250 meter buffer -12.11000  0.09448     -0.30210  -0.7732   9.0200 15 -176.746
## 4500 meter buffer -20.69000  0.58330      0.50590  -3.3870  22.5300 15 -178.915
## 4250 meter buffer -17.77000 -0.08114      0.09835  -3.6280  17.9600 15 -179.394
## 3000 meter buffer  -9.50600  0.35300     -0.08657  -0.2605   8.4640 15 -179.732
## 3500 meter buffer -12.67000 -1.09100     -1.27100  -1.8420  11.0800 15 -180.010
## 2750 meter buffer  -9.04700 -0.05575     -0.36520  -0.3487   7.4750 15 -181.317
## 4750 meter buffer -17.72000  2.34000      1.89700  -1.2780  21.2300 15 -184.111
## 5000 meter buffer -15.40000  2.43200      1.84400  -0.3963  19.0600 15 -185.542
## 2500 meter buffer  -7.58600  0.58260      0.11670   0.6559   4.7480 15 -186.923
## 1750 meter buffer  -4.38400  3.00900      2.39100   1.5010   2.4770 15 -187.397
## 1250 meter buffer  -3.15700  3.50500      2.39800   1.4510   1.2450 15 -187.748
## 1500 meter buffer  -3.96700  3.33900      2.32400   1.5970   3.0150 15 -188.914
## 2000 meter buffer  -5.03600  2.70000      1.97100   1.1880   4.2230 15 -190.436
## 2250 meter buffer  -6.08400  1.28100      0.73260   0.8988   2.7970 15 -190.880
## 1000 meter buffer   0.03988  0.36850      1.22600   1.5110  -1.8270 15 -194.098
## 750 meter buffer    0.43300  1.77100      1.41300   1.2110  -0.6427 15 -200.541
## 500 meter buffer    0.80290  1.75000      0.98560   1.0630  -0.7946 15 -200.707
##                    AICc delta weight
## 3750 meter buffer 383.1  0.00  0.678
## 4000 meter buffer 385.6  2.47  0.197
## 3250 meter buffer 387.0  3.89  0.097
## 4500 meter buffer 391.4  8.22  0.011
## 4250 meter buffer 392.3  9.18  0.007
## 3000 meter buffer 393.0  9.86  0.005
## 3500 meter buffer 393.5 10.41  0.004
## 2750 meter buffer 396.2 13.03  0.001
## 4750 meter buffer 401.8 18.62  0.000
## 5000 meter buffer 404.6 21.48  0.000
## 2500 meter buffer 407.4 24.24  0.000
## 1750 meter buffer 408.3 25.19  0.000
## 1250 meter buffer 409.0 25.89  0.000
## 1500 meter buffer 411.4 28.22  0.000
## 2000 meter buffer 414.4 31.27  0.000
## 2250 meter buffer 415.3 32.16  0.000
## 1000 meter buffer 421.7 38.59  0.000
## 750 meter buffer  434.6 51.48  0.000
## 500 meter buffer  434.9 51.81  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For grey wolf top model is 4500 m buffer

Let’s get the model summary for this model

summary(wolf_mods_no250$`4500 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(grey_wolf, absent_grey_wolf) ~ seismic_lines + pipeline +  
##     borrowpits + wellsites + roads + trails + lc_class20 + lc_class34 +  
##     lc_class50 + lc_class110 + lc_class210 + lc_class220 + lc_class230 +  
##     (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    387.8    433.2   -178.9    357.8      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  array  (Intercept) 4.188e-09 6.472e-05
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -6.9362     5.0762  -1.366   0.1718    
## seismic_lines   0.5059     1.6039   0.315   0.7525    
## pipeline      -20.6864     3.8734  -5.341 9.26e-08 ***
## borrowpits     12.2580     9.7501   1.257   0.2087    
## wellsites      22.5258     4.0771   5.525 3.29e-08 ***
## roads           0.5833     1.9442   0.300   0.7641    
## trails         -3.3871     1.9879  -1.704   0.0884 .  
## lc_class20      9.3963     6.4449   1.458   0.1449    
## lc_class34    -83.9732    26.9038  -3.121   0.0018 ** 
## lc_class50      4.2511     4.8758   0.872   0.3833    
## lc_class110     5.4079     5.0393   1.073   0.2832    
## lc_class210     8.5561     5.1606   1.658   0.0973 .  
## lc_class220     2.4732     4.9505   0.500   0.6174    
## lc_class230     3.4749     4.9971   0.695   0.4868    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lc_class34 still presenting some issues, interesting that seismic lines weren’t significant and have a negative estimate

Lynx

lynx_mods_no250 <- final_df %>%
  
  # remove 350 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(lynx, absent_lynx) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(lynx_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 750 meter buffer     -7.2440          +   4.1720      5.3780      -3.860
## 1000 meter buffer   -10.3500          +   2.0140      8.1860       2.235
## 1250 meter buffer   -13.7800          +  -1.9290     10.7600      12.920
## 500 meter buffer     -1.6740          +  -4.5310     -0.1358      -4.147
## 2500 meter buffer   -15.9300          +  -9.7760     12.2600      14.020
## 1750 meter buffer   -14.3300          +   1.0280     11.1100      14.550
## 1500 meter buffer   -15.0400          +  -0.3430     11.3900      12.660
## 5000 meter buffer    -2.3840          +  -0.6378     -0.8068      -0.644
## 2750 meter buffer   -14.3000          +  -9.0810     10.9000      12.640
## 4750 meter buffer    -0.1961          +   0.5363     -2.8660      -4.934
## 2000 meter buffer   -14.0100          +  -0.5195     10.6900      11.260
## 3250 meter buffer    -6.6590          +  -4.9000      4.0480       3.626
## 2250 meter buffer   -14.1600          +  -4.7400     11.2800      10.870
## 3500 meter buffer    -4.5340          +  -0.8624      1.8270       1.040
## 3000 meter buffer    -9.9420          +  -5.7680      7.1950       7.039
## 4500 meter buffer    -2.1510          +  -0.2477     -0.9033      -2.676
## 3750 meter buffer    -2.9910          +  -1.2510      0.6628      -1.126
## 4000 meter buffer    -0.6971          +  -1.1460     -0.6745      -5.145
## 4250 meter buffer    -0.5769          +  -0.6674     -1.3080      -5.281
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 750 meter buffer       5.7510      4.8420      3.5120      3.6320      3.7990
## 1000 meter buffer      8.1460      7.4930      6.4920      8.3480      6.8070
## 1250 meter buffer     11.5600     10.2800      9.1820      9.0880      9.0930
## 500 meter buffer       0.5488     -0.6783     -1.1560     -1.5890     -1.3020
## 2500 meter buffer     11.8900     10.8500      9.6450      2.2140     10.3800
## 1750 meter buffer     11.8000     10.3500      9.4540      4.1960      8.9840
## 1500 meter buffer     13.4700     11.4900      9.6920      9.5520     10.4900
## 5000 meter buffer     -0.1329     -1.1260     -3.8710     25.3300     -1.5540
## 2750 meter buffer     10.3800     10.3500      9.1990      8.3780     10.4100
## 4750 meter buffer     -2.4280     -2.8490     -5.3280     20.5100     -3.4740
## 2000 meter buffer     10.9300      9.1630      8.8460      3.1450      8.0030
## 3250 meter buffer      2.5390      3.4900      1.8720     11.9400      3.8980
## 2250 meter buffer     10.4000      9.6760      9.3810      0.3421      8.9210
## 3500 meter buffer      1.2710      1.4870     -0.3037     11.6000      1.6080
## 3000 meter buffer      5.2430      6.9300      4.9650     10.1700      6.6630
## 4500 meter buffer     -0.3603     -1.5210     -2.9030     18.0600     -1.0670
## 3750 meter buffer     -0.2613     -0.3054     -1.4640      9.3620      0.1612
## 4000 meter buffer     -2.6760     -2.2500     -2.3350     11.4500     -1.8060
## 4250 meter buffer     -2.2900     -3.0090     -3.2680     11.1100     -2.0510
##                   cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df   logLik
## 750 meter buffer    1.0810   2.0560      0.48590   0.4565   1.1230 15 -288.714
## 1000 meter buffer   2.5290   2.2420      0.91990   1.1600  -1.0590 15 -293.455
## 1250 meter buffer   3.3230   2.8570      1.61800   2.1720  -0.9282 15 -293.460
## 500 meter buffer   -0.6344   0.9736      0.04477  -0.1048   2.9010 15 -295.259
## 2500 meter buffer   3.2780   4.5130      3.82400   2.9140   3.8330 15 -295.629
## 1750 meter buffer   2.3060   3.0990      2.63300   2.3020   1.2800 15 -295.751
## 1500 meter buffer   2.2210   2.9480      1.93200   1.7180   2.0870 15 -296.296
## 5000 meter buffer  -6.6520   2.5920      2.29900  -0.7819  17.4400 15 -296.903
## 2750 meter buffer   3.0290   3.0410      2.57200   1.0790   5.0200 15 -297.569
## 4750 meter buffer  -6.1050   2.1790      1.97500  -0.6911  15.1400 15 -298.569
## 2000 meter buffer   2.2040   3.8400      3.30900   3.3810   1.9170 15 -298.671
## 3250 meter buffer  -1.4270   2.0360      2.12900  -0.7995   8.9050 15 -299.136
## 2250 meter buffer   3.9250   3.5440      2.83800   3.0790   1.8860 15 -299.400
## 3500 meter buffer  -2.2590   1.7650      2.04700  -1.1090   9.6540 15 -301.128
## 3000 meter buffer   2.0190   2.3400      2.26800   0.5232   4.8820 15 -301.414
## 4500 meter buffer  -3.3690   2.0510      1.96500  -0.5899  11.2200 15 -302.073
## 3750 meter buffer  -2.2100   1.6070      1.86400  -0.9372   9.3990 15 -302.805
## 4000 meter buffer  -1.2260   0.6927      1.14300  -1.2400   7.7790 15 -302.955
## 4250 meter buffer  -1.3800   1.1680      1.56700  -1.0140   7.7080 15 -303.274
##                    AICc delta weight
## 750 meter buffer  611.0  0.00  0.979
## 1000 meter buffer 620.4  9.48  0.009
## 1250 meter buffer 620.4  9.49  0.009
## 500 meter buffer  624.0 13.09  0.001
## 2500 meter buffer 624.8 13.83  0.001
## 1750 meter buffer 625.0 14.07  0.001
## 1500 meter buffer 626.1 15.16  0.000
## 5000 meter buffer 627.3 16.38  0.000
## 2750 meter buffer 628.7 17.71  0.000
## 4750 meter buffer 630.7 19.71  0.000
## 2000 meter buffer 630.9 19.91  0.000
## 3250 meter buffer 631.8 20.84  0.000
## 2250 meter buffer 632.3 21.37  0.000
## 3500 meter buffer 635.8 24.83  0.000
## 3000 meter buffer 636.4 25.40  0.000
## 4500 meter buffer 637.7 26.72  0.000
## 3750 meter buffer 639.1 28.18  0.000
## 4000 meter buffer 639.4 28.48  0.000
## 4250 meter buffer 640.1 29.12  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For lynx the top model is the 1000 m buffer

Let’s get the model summary

summary(lynx_mods_no250$`1000 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(lynx, absent_lynx) ~ seismic_lines + pipeline + borrowpits +  
##     wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +  
##     lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 |      array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    616.9    662.3   -293.5    586.9      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.009742 0.0987  
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -10.3480     2.6475  -3.909 9.28e-05 ***
## seismic_lines   0.9199     0.5963   1.543  0.12289    
## pipeline        2.5293     0.9589   2.638  0.00835 ** 
## borrowpits      2.0140     4.1551   0.485  0.62789    
## wellsites      -1.0595     1.3301  -0.797  0.42571    
## roads           2.2419     0.8048   2.786  0.00534 ** 
## trails          1.1598     0.6942   1.671  0.09478 .  
## lc_class20      2.2347     3.5463   0.630  0.52860    
## lc_class34      8.3482     3.1569   2.644  0.00818 ** 
## lc_class50      6.8067     2.6463   2.572  0.01011 *  
## lc_class110     8.1856     2.5943   3.155  0.00160 ** 
## lc_class210     8.1459     2.7213   2.993  0.00276 ** 
## lc_class220     7.4935     2.6936   2.782  0.00540 ** 
## lc_class230     6.4919     2.6612   2.439  0.01471 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Moose

moose_mods_no250 <- final_df %>%
  
  # remove 350 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(moose, absent_moose) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(moose_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 3750 meter buffer   -0.18680          +   1.6060     -1.6770     -1.5320
## 500 meter buffer    -5.31100          +  -8.7360      4.9570      5.4540
## 4000 meter buffer   -0.06881          +  -3.5510     -1.6430     -2.0670
## 3250 meter buffer   -1.86900          +   8.3310      0.7309     -0.9436
## 3000 meter buffer   -0.85170          +   9.3060     -0.1209     -0.3379
## 2250 meter buffer   -0.84560          +  10.1400      1.0920      0.8478
## 3500 meter buffer   -0.08749          +   3.5810     -1.4530     -1.1310
## 750 meter buffer    -4.78000          +  -0.7357      2.9290      7.0660
## 2750 meter buffer   -0.17090          +  11.7400      0.1404     -1.0460
## 2500 meter buffer   -0.18740          +  11.5200      0.3941     -0.4472
## 4750 meter buffer    0.88060          + -10.3800     -2.7910     -3.6680
## 5000 meter buffer    2.04700          +  -9.2220     -3.4390     -4.4980
## 4500 meter buffer    0.58590          +  -6.9310     -2.2090     -3.4300
## 4250 meter buffer    0.89340          +  -4.4540     -2.3180     -4.1220
## 2000 meter buffer   -0.90940          +   4.3410      0.1203      0.5869
## 1750 meter buffer   -1.74400          +   3.7490      1.2330      1.9700
## 1500 meter buffer   -0.04846          +  -1.8850     -0.3197      0.9663
## 1250 meter buffer   -0.97140          +  -3.1900      0.3438      1.9300
## 1000 meter buffer   -1.12700          +  -1.7150      0.4704      2.3250
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 3750 meter buffer    -0.74180     -1.4910      0.1646    -19.5100    -1.40100
## 500 meter buffer      4.00100      5.6370      5.6360      4.5460     4.02900
## 4000 meter buffer     0.05318     -2.1260     -0.2459     -1.6980    -1.76700
## 3250 meter buffer     1.44200      0.7038      2.0560    -12.4800     0.21500
## 3000 meter buffer     0.43970     -0.7143      0.7540    -15.2600    -0.80070
## 2250 meter buffer     0.75530     -1.3000      1.5580    -15.6500    -0.42090
## 3500 meter buffer    -0.70580     -1.3730      0.0586    -18.1800    -1.53500
## 750 meter buffer      4.06000      4.5610      4.1730      2.3200     2.94100
## 2750 meter buffer     0.05767     -1.3620      0.3643    -19.2800    -1.30400
## 2500 meter buffer     0.18610     -1.9870      0.1478    -18.5900    -1.43000
## 4750 meter buffer    -0.62150     -3.0420     -0.4440     18.9700    -2.66500
## 5000 meter buffer    -1.57000     -3.9060     -1.4270     21.9600    -3.53000
## 4500 meter buffer    -0.77680     -2.6810     -0.3695      4.7680    -2.45700
## 4250 meter buffer    -0.81380     -2.9140     -0.5764      1.7190    -2.32900
## 2000 meter buffer     0.03522     -0.7530      1.4490     -8.1510    -0.70010
## 1750 meter buffer     0.96740      0.4748      1.9070     -4.1270    -0.22920
## 1500 meter buffer    -0.79890     -0.2145     -0.5274     -2.9790    -1.56800
## 1250 meter buffer    -0.19010      0.6056      0.6409     -0.6084    -0.43920
## 1000 meter buffer    -0.24210      1.0790      1.0800      0.6271    -0.01273
##                   cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df   logLik
## 3750 meter buffer   -6.291  -0.5856      1.00500 -0.09781   0.5362 15 -302.076
## 500 meter buffer    -1.439  -1.2940     -0.64690 -0.98350  -0.4259 15 -302.597
## 4000 meter buffer   -6.063  -0.2953      1.03600 -0.01724  -1.0620 15 -302.787
## 3250 meter buffer   -6.637  -1.1550      0.43790 -0.32210   1.0850 15 -303.233
## 3000 meter buffer   -7.115  -1.2020      0.60440 -0.82810   3.1200 15 -303.251
## 2250 meter buffer   -5.015  -1.7170     -0.01713 -1.09600   1.1260 15 -303.374
## 3500 meter buffer   -7.034  -0.7724      0.81050 -0.49360   2.0380 15 -303.837
## 750 meter buffer    -1.838  -0.9000     -0.01174 -1.10100  -0.8452 15 -303.986
## 2750 meter buffer   -6.737  -1.5230      0.16060 -1.12800   3.0270 15 -303.995
## 2500 meter buffer   -5.769  -1.3710      0.30700 -0.75810   1.7990 15 -304.601
## 4750 meter buffer   -5.573  -0.5700      0.66620 -0.08074  -0.5954 15 -304.864
## 5000 meter buffer   -5.110  -1.2730      0.41010 -0.57930  -1.4990 15 -304.984
## 4500 meter buffer   -4.892  -0.4925      0.87600  0.19800  -2.0140 15 -306.058
## 4250 meter buffer   -4.950  -0.8266      0.55860 -0.13530  -2.2470 15 -306.128
## 2000 meter buffer   -4.541  -1.1100      0.29120 -0.65190   0.7489 15 -306.878
## 1750 meter buffer   -3.957  -1.0500      0.26310 -0.35570  -1.1090 15 -307.198
## 1500 meter buffer   -4.002  -1.2080     -0.02248 -0.91600  -0.5908 15 -308.780
## 1250 meter buffer   -2.796  -1.1330     -0.11170 -0.84400  -1.4450 15 -311.149
## 1000 meter buffer   -2.204  -1.5150     -0.38660 -1.03500  -1.0020 15 -313.019
##                    AICc delta weight
## 3750 meter buffer 637.7  0.00  0.271
## 500 meter buffer  638.7  1.04  0.161
## 4000 meter buffer 639.1  1.42  0.133
## 3250 meter buffer 640.0  2.31  0.085
## 3000 meter buffer 640.0  2.35  0.084
## 2250 meter buffer 640.3  2.60  0.074
## 3500 meter buffer 641.2  3.52  0.047
## 750 meter buffer  641.5  3.82  0.040
## 2750 meter buffer 641.5  3.84  0.040
## 2500 meter buffer 642.7  5.05  0.022
## 4750 meter buffer 643.3  5.58  0.017
## 5000 meter buffer 643.5  5.82  0.015
## 4500 meter buffer 645.6  7.97  0.005
## 4250 meter buffer 645.8  8.10  0.005
## 2000 meter buffer 647.3  9.60  0.002
## 1750 meter buffer 647.9 10.24  0.002
## 1500 meter buffer 651.1 13.41  0.000
## 1250 meter buffer 655.8 18.15  0.000
## 1000 meter buffer 659.6 21.89  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For moose the top model is the 3750 m buffer

Let’s get the model summary

summary(moose_mods_no250$`3750 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(moose, absent_moose) ~ seismic_lines + pipeline + borrowpits +  
##     wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +  
##     lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 |      array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    634.2    679.5   -302.1    604.2      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev.
##  array  (Intercept) 4.409e-10 2.1e-05 
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -0.18680    2.73864  -0.068 0.945619    
## seismic_lines   1.00450    0.88172   1.139 0.254598    
## pipeline       -6.29073    1.65572  -3.799 0.000145 ***
## borrowpits      1.60628    5.40801   0.297 0.766452    
## wellsites       0.53622    2.16428   0.248 0.804319    
## roads          -0.58556    1.05508  -0.555 0.578901    
## trails         -0.09781    0.99686  -0.098 0.921839    
## lc_class20     -1.53153    3.35843  -0.456 0.648373    
## lc_class34    -19.51133   10.62742  -1.836 0.066366 .  
## lc_class50     -1.40080    2.61503  -0.536 0.592186    
## lc_class110    -1.67652    2.64796  -0.633 0.526646    
## lc_class210    -0.74176    2.85993  -0.259 0.795354    
## lc_class220    -1.49130    2.66783  -0.559 0.576168    
## lc_class230     0.16455    2.59309   0.063 0.949401    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Red fox

fox_mods_no250 <- final_df %>%
  
  # remove 350 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(red_fox, absent_red_fox) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(fox_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 500 meter buffer     -2.5350          +    8.045     0.30030     -4.2490
## 3250 meter buffer    -1.8780          +  -24.310     6.31300      0.3432
## 5000 meter buffer    -4.2450          +  -48.150    -0.01949     15.2500
## 4750 meter buffer    -6.6770          +  -50.640     2.17400     17.0500
## 3000 meter buffer    -3.2410          +  -36.230     8.11300     -0.6053
## 3500 meter buffer    -1.0710          +  -28.270     6.55400     -2.3710
## 4500 meter buffer    -5.0310          +  -45.170     3.16900     13.1900
## 4250 meter buffer    -5.3540          +  -36.060     3.83200     13.2200
## 3750 meter buffer    -5.3600          +  -44.530     5.79400      6.1910
## 4000 meter buffer    -6.7640          +  -34.860     5.17100      7.6620
## 2750 meter buffer    -4.4030          +  -34.330     6.49000      1.7550
## 750 meter buffer     -4.8760          +   -7.385     1.66500      0.2266
## 1000 meter buffer    -0.3164          +    1.064    -3.20400     -3.8810
## 2250 meter buffer    -7.0970          +  -35.320     4.10400      2.7390
## 2500 meter buffer    -2.7410          +  -20.620     3.61600     -0.1645
## 2000 meter buffer    -6.0340          +  -20.030     4.07300      2.7720
## 1250 meter buffer    -1.7950          +   -8.857    -2.70400     -4.4450
## 1500 meter buffer    -3.4350          +   -5.641     0.97020     -1.8670
## 1750 meter buffer    -3.3080          +  -14.860     0.86500     -2.2250
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 500 meter buffer      -7.5030     -3.8590     -4.6400     -5.5180     -5.2100
## 3250 meter buffer      5.7970     -3.8260     -0.1972      1.3640     -2.3330
## 5000 meter buffer      5.9290     -5.3100     -6.6610      5.8790     -4.2780
## 4750 meter buffer      6.5070     -3.3150     -4.0770     31.0900     -1.7330
## 3000 meter buffer      8.3460     -1.7820      0.4510     -7.3680     -1.7250
## 3500 meter buffer      3.9730     -4.1260      0.3306      9.0720     -2.3890
## 4500 meter buffer      5.5000     -3.3350     -4.6490     10.2500     -2.9750
## 4250 meter buffer      5.9440     -3.4510     -2.9620     23.9100     -2.3760
## 3750 meter buffer      4.9450     -1.6920     -0.3727     28.0600     -0.8399
## 4000 meter buffer      5.8850     -1.4440     -2.4020     32.9100     -0.5652
## 2750 meter buffer      7.6970     -0.7962      4.5920     -0.8732      1.0240
## 750 meter buffer      -3.5300      1.0090      0.2439     -0.9965      1.1890
## 1000 meter buffer     -6.3970     -6.8690     -5.2280     -4.5860     -4.3310
## 2250 meter buffer      4.9280      3.3350      4.1270     16.9000      2.7110
## 2500 meter buffer      2.9930     -1.1920      1.6660     -2.1310     -0.5620
## 2000 meter buffer      0.1995      2.8610      3.3670     11.0100      2.6440
## 1250 meter buffer     -4.0600     -4.2180     -3.4020      1.4780     -2.8780
## 1500 meter buffer     -2.0290     -1.1090     -0.9122      3.3400      0.2728
## 1750 meter buffer     -0.5451     -1.0060     -1.8550      5.4800      0.1459
##                   cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df   logLik
## 500 meter buffer     5.077  5.75000       3.9700  1.92000 -14.8600 15 -113.738
## 3250 meter buffer    7.644 -3.02600      -2.7870 -1.96100 -12.7800 15 -134.069
## 5000 meter buffer    4.249  4.23400       5.7080 -2.09700   0.6589 15 -134.197
## 4750 meter buffer    7.789  4.74200       5.8910 -0.41150  -5.8730 15 -134.722
## 3000 meter buffer    3.346 -2.58400      -2.9190 -2.14600  -4.0520 15 -134.936
## 3500 meter buffer    8.457 -3.99700      -3.1680 -2.06200 -15.9000 15 -135.792
## 4500 meter buffer    4.204  3.40000       3.9740 -0.71710  -1.5480 15 -136.499
## 4250 meter buffer    7.043  2.81500       2.8250  1.02400  -9.1310 15 -136.936
## 3750 meter buffer    8.756  1.02400       1.0380  1.21400 -12.9400 15 -137.141
## 4000 meter buffer    9.758  3.37400       3.1730  2.54600 -15.8300 15 -138.440
## 2750 meter buffer    4.948 -4.11700      -3.7050 -2.44900  -3.1040 15 -138.678
## 750 meter buffer     2.669  2.89800       1.3630  0.17270   3.8360 15 -139.781
## 1000 meter buffer    1.756  3.24100       1.8410 -1.06900   3.8370 15 -142.960
## 2250 meter buffer    6.604 -0.59570      -1.2960 -0.03987  -2.8970 15 -145.048
## 2500 meter buffer    5.403 -3.45200      -2.9160 -1.96200  -3.6040 15 -145.241
## 2000 meter buffer    8.459 -1.19700      -1.4540  1.12200  -8.6090 15 -145.514
## 1250 meter buffer    1.103  3.06200       0.6782  1.82400   2.2090 15 -150.979
## 1500 meter buffer    2.857  0.21900      -1.0510  1.86200  -2.1440 15 -152.307
## 1750 meter buffer    1.678 -0.01589      -0.8923  0.40610   1.0990 15 -153.758
##                    AICc delta weight
## 500 meter buffer  261.0  0.00      1
## 3250 meter buffer 301.7 40.66      0
## 5000 meter buffer 301.9 40.92      0
## 4750 meter buffer 303.0 41.97      0
## 3000 meter buffer 303.4 42.40      0
## 3500 meter buffer 305.1 44.11      0
## 4500 meter buffer 306.5 45.52      0
## 4250 meter buffer 307.4 46.40      0
## 3750 meter buffer 307.8 46.81      0
## 4000 meter buffer 310.4 49.40      0
## 2750 meter buffer 310.9 49.88      0
## 750 meter buffer  313.1 52.09      0
## 1000 meter buffer 319.5 58.45      0
## 2250 meter buffer 323.6 62.62      0
## 2500 meter buffer 324.0 63.01      0
## 2000 meter buffer 324.6 63.55      0
## 1250 meter buffer 335.5 74.48      0
## 1500 meter buffer 338.1 77.14      0
## 1750 meter buffer 341.0 80.04      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For red fox the top model is 3750 m buffer

Let’s get the model summary

summary(fox_mods_no250$`3750 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(red_fox, absent_red_fox) ~ seismic_lines + pipeline + borrowpits +  
##     wellsites + roads + trails + lc_class20 + lc_class34 + lc_class50 +  
##     lc_class110 + lc_class210 + lc_class220 + lc_class230 + (1 |      array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    304.3    349.6   -137.1    274.3      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  array  (Intercept) 7.606e-10 2.758e-05
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)    -5.3604     5.6689  -0.946  0.34437   
## seismic_lines   1.0379     2.4927   0.416  0.67715   
## pipeline        8.7560     3.4361   2.548  0.01083 * 
## borrowpits    -44.5280    16.6486  -2.675  0.00748 **
## wellsites     -12.9410     6.4374  -2.010  0.04440 * 
## roads           1.0242     2.9737   0.344  0.73052   
## trails          1.2141     2.8278   0.429  0.66767   
## lc_class20      6.1909     6.7804   0.913  0.36122   
## lc_class34     28.0557    21.5870   1.300  0.19372   
## lc_class50     -0.8399     5.1401  -0.163  0.87020   
## lc_class110     5.7942     5.0389   1.150  0.25018   
## lc_class210     4.9452     5.7475   0.860  0.38956   
## lc_class220    -1.6918     5.3634  -0.315  0.75243   
## lc_class230    -0.3727     4.9820  -0.075  0.94037   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Gagh! Borrow pits does not have a reasonable estimate and SE

White-tailed deer

deer_mods_no250 <- final_df %>%
  
  # remove 350 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

      # have to include the `` around the white-tailed_deer or R won't recognize it as a variable because of the -
     glmmTMB::glmmTMB(cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~
                        seismic_lines +
                        pipeline +
                        borrowpits +
                        wellsites +
                        roads +
                        trails + 
                        lc_class20 +
                        lc_class34 +
                        lc_class50 +
                        lc_class110 +
                        lc_class210 +
                        lc_class220 +
                        lc_class230 +
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(deer_mods_no250)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(brr) cnd(lc_c11) cnd(lc_c20)
## 2500 meter buffer      4.449          + -13.9700     -2.9190     -1.7080
## 2750 meter buffer      6.563          + -12.9600     -4.4960     -3.1680
## 4750 meter buffer      2.988          + -11.1400     -0.2994     -6.3250
## 5000 meter buffer      3.271          + -10.1700     -0.2263     -7.7080
## 2250 meter buffer      3.505          + -12.0900     -2.2360     -0.8935
## 4500 meter buffer      1.816          + -11.2500      0.4433     -4.8170
## 4250 meter buffer      2.148          + -15.2500     -0.5879     -3.2000
## 3000 meter buffer      3.776          + -16.1400     -2.6610     -2.0220
## 3750 meter buffer      2.684          + -15.4900     -1.2780     -3.4710
## 4000 meter buffer      1.762          + -18.9700     -0.8749     -2.1220
## 2000 meter buffer      3.126          + -10.4100     -2.0760     -1.5350
## 3250 meter buffer      4.144          +  -5.8400     -2.2140     -3.8770
## 3500 meter buffer      5.690          +  -6.3760     -3.2390     -5.3640
## 1750 meter buffer      3.550          +  -9.0480     -2.4720     -2.0090
## 1500 meter buffer      3.826          + -13.8900     -2.5630     -3.4080
## 750 meter buffer       1.911          +  -5.3960     -1.5810     -1.1940
## 1250 meter buffer      4.386          +  -8.2400     -2.9200     -3.6650
## 1000 meter buffer      4.165          + -10.6500     -2.7420     -3.0780
## 500 meter buffer       1.976          +  -0.3651     -0.7439     -0.8363
##                   cnd(lc_c21) cnd(lc_c22) cnd(lc_c23) cnd(lc_c34) cnd(lc_c50)
## 2500 meter buffer    -0.87940    -1.46300     -4.2760     -59.850      -5.517
## 2750 meter buffer    -2.37500    -3.51700     -5.9550     -71.510      -7.653
## 4750 meter buffer     2.00200    -0.06123     -0.9652    -149.000      -5.193
## 5000 meter buffer     1.23400     0.22020     -0.7564    -165.500      -5.356
## 2250 meter buffer     0.09957     0.14690     -3.5810     -43.730      -3.789
## 4500 meter buffer     3.47500    -0.13640     -0.7182    -131.100      -4.043
## 4250 meter buffer     2.39300    -0.33740     -1.7970    -120.000      -4.636
## 3000 meter buffer    -0.16560    -0.91040     -4.9410     -66.660      -5.607
## 3750 meter buffer     2.13000    -1.15300     -3.8980     -81.490      -5.012
## 4000 meter buffer     3.29000    -0.66540     -2.9250     -95.770      -4.413
## 2000 meter buffer    -0.31490    -0.17840     -3.2490     -32.570      -3.927
## 3250 meter buffer    -0.09955    -1.34500     -4.3950     -77.180      -5.775
## 3500 meter buffer    -0.93680    -2.97200     -5.2020     -88.770      -7.019
## 1750 meter buffer    -1.77900    -1.19900     -2.8100     -25.590      -4.363
## 1500 meter buffer    -3.39300    -1.80200     -3.4400     -18.810      -5.015
## 750 meter buffer     -0.87950     0.86670     -0.5545      -6.163      -1.193
## 1250 meter buffer    -5.14100    -2.82300     -2.9750     -16.550      -4.948
## 1000 meter buffer    -4.75600    -2.13700     -2.2190     -12.210      -3.631
## 500 meter buffer     -1.13800     0.36320     -0.1688      -4.073      -1.055
##                   cnd(ppl) cnd(rds) cnd(ssm_lns) cnd(trl) cnd(wll) df   logLik
## 2500 meter buffer  0.48160 -1.32300      -1.2990 -2.08000   2.8420 15 -336.320
## 2750 meter buffer -0.58130 -1.78500      -1.6890 -2.29900   4.1070 15 -337.058
## 4750 meter buffer  0.86930 -1.52400      -1.9100  0.03205  -1.9650 15 -338.898
## 5000 meter buffer  0.46840 -1.68600      -2.1540 -0.10530  -2.0670 15 -339.035
## 2250 meter buffer  0.55350 -1.62300      -1.8310 -2.53700   2.1940 15 -340.144
## 4500 meter buffer  1.02900 -0.92010      -1.2660  0.27020  -1.9190 15 -340.568
## 4250 meter buffer  0.54160 -0.42630      -0.7360  0.24770  -0.1690 15 -341.224
## 3000 meter buffer -0.93460 -0.37020      -0.7497 -1.23900   4.4660 15 -341.911
## 3750 meter buffer -0.63730 -0.17760      -0.5014 -0.40380   2.5760 15 -344.834
## 4000 meter buffer  0.02576  0.34870      -0.1156  0.27040   0.5567 15 -345.197
## 2000 meter buffer  0.87530 -1.14000      -1.4120 -1.89000   1.9070 15 -345.887
## 3250 meter buffer -1.22900 -1.11800      -1.1650 -1.32200   2.5980 15 -346.669
## 3500 meter buffer -1.91000 -1.87000      -1.6590 -1.91900   3.7660 15 -347.851
## 1750 meter buffer  0.43570 -0.85500      -1.1890 -1.75800   0.8253 15 -353.649
## 1500 meter buffer  0.60190 -0.01901      -0.7268 -0.57980  -0.9992 15 -360.550
## 750 meter buffer   0.23730 -1.47500      -1.8050 -2.04200   0.1912 15 -367.854
## 1250 meter buffer  0.31870 -0.21480      -0.5859 -0.63250  -1.6450 15 -368.323
## 1000 meter buffer  0.05693 -0.89100      -1.0910 -1.37100  -1.4910 15 -370.040
## 500 meter buffer   0.45850 -2.20500      -2.0600 -1.11900  -0.9517 15 -373.158
##                    AICc delta weight
## 2500 meter buffer 706.2  0.00  0.600
## 2750 meter buffer 707.6  1.48  0.287
## 4750 meter buffer 711.3  5.16  0.045
## 5000 meter buffer 711.6  5.43  0.040
## 2250 meter buffer 713.8  7.65  0.013
## 4500 meter buffer 714.7  8.50  0.009
## 4250 meter buffer 716.0  9.81  0.004
## 3000 meter buffer 717.4 11.18  0.002
## 3750 meter buffer 723.2 17.03  0.000
## 4000 meter buffer 723.9 17.76  0.000
## 2000 meter buffer 725.3 19.14  0.000
## 3250 meter buffer 726.9 20.70  0.000
## 3500 meter buffer 729.2 23.06  0.000
## 1750 meter buffer 740.8 34.66  0.000
## 1500 meter buffer 754.6 48.46  0.000
## 750 meter buffer  769.2 63.07  0.000
## 1250 meter buffer 770.2 64.01  0.000
## 1000 meter buffer 773.6 67.44  0.000
## 500 meter buffer  779.8 73.68  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For deer the top model was also the 3750 buffer

Let’s get the model summary

summary(deer_mods_no250$`3750 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~ seismic_lines +  
##     pipeline + borrowpits + wellsites + roads + trails + lc_class20 +  
##     lc_class34 + lc_class50 + lc_class110 + lc_class210 + lc_class220 +  
##     lc_class230 + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    719.7    765.0   -344.8    689.7      137 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 1.021    1.01    
## Number of obs: 152, groups:  array, 4
## 
## Conditional model:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     2.6845     2.8508   0.942  0.34636    
## seismic_lines  -0.5014     0.9714  -0.516  0.60574    
## pipeline       -0.6373     1.6032  -0.398  0.69096    
## borrowpits    -15.4939     5.8852  -2.633  0.00847 ** 
## wellsites       2.5761     2.5555   1.008  0.31342    
## roads          -0.1776     1.2096  -0.147  0.88329    
## trails         -0.4038     1.1031  -0.366  0.71434    
## lc_class20     -3.4711     3.4728  -1.000  0.31755    
## lc_class34    -81.4906    15.1035  -5.395 6.83e-08 ***
## lc_class50     -5.0125     2.6098  -1.921  0.05478 .  
## lc_class110    -1.2782     2.7120  -0.471  0.63742    
## lc_class210     2.1301     2.8871   0.738  0.46064    
## lc_class220    -1.1534     2.6142  -0.441  0.65907    
## lc_class230    -3.8981     2.6216  -1.487  0.13703    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Coefficient plots

This is also a dog shit way to do this but I need to get this done

coeffs <- confint(black_bear_mods_no250$`500 meter buffer`) %>% 
  
  as_tibble() %>%
  
  # subset to just the HFI variables for these plots
    slice_head(n = 6) %>% 
  
  # add a column where we can put the feature names
  rowid_to_column() %>% 
  
  # rename the columns
  rename('lower' = `2.5 %`,
         'upper' = `97.5 %`,
         'estimate' = Estimate,
         'feature' = rowid) %>% 
  
  # rename the entries to features, need to look at the order the features are in from the model summary and ensure it matches
  mutate(feature = as.factor(feature),
         feature = recode(feature,
                          '1' = 'seismic_lines',
                          '2' = 'pipeline',
                          '3' = 'borrowpits',
                          '4' = 'wellsites',
                          '5' = 'roads',
                          '6' = 'trails'))
  
# print
 coeffs
## # A tibble: 6 × 4
##   feature        lower upper estimate
##   <fct>          <dbl> <dbl>    <dbl>
## 1 seismic_lines -3.52  2.78    -0.370
## 2 pipeline      -0.863 0.615   -0.124
## 3 borrowpits    -1.53  0.731   -0.401
## 4 wellsites     -9.74  4.28    -2.73 
## 5 roads         -1.12  1.98     0.430
## 6 trails        -1.14  0.883   -0.131

Try ad do this for all top models at once with purrr

top_models <- list(black_bear_mods_no250$`500 meter buffer`,
                   caribou_mods_no250$`1250 meter buffer`,
                   coyote_mods_no250$`4500 meter buffer`,
                   fisher_mods_no250$`2000 meter buffer`,
                   wolf_mods_no250$`4500 meter buffer`,
                   lynx_mods_no250$`1000 meter buffer`,
                   moose_mods_no250$`3750 meter buffer`,
                   fox_mods_no250$`3750 meter buffer`,
                   deer_mods_no250$`3750 meter buffer`) %>%  
  
  # use purrr to create coefficient table for all models
  purrr::map(
    ~.x %>% 
      
      confint() %>% 
 
      as_tibble() %>%
      
      # subset to just the HFI variables for these plots
      slice_head(n = 6) %>% 
      
      # add a column where we can put the feature names
      rowid_to_column() %>% 
      
      # rename the columns
      rename('lower' = `2.5 %`,
             'upper' = `97.5 %`,
             'estimate' = Estimate,
             'feature' = rowid) %>% 
      
      # rename the entries to features, need to look at the order the features are in from the model summary and ensure it matches
      mutate(feature = as.factor(feature),
             feature = recode(feature,
                              '1' = 'seismic_lines',
                              '2' = 'pipeline',
                              '3' = 'borrowpits',
                              '4' = 'wellsites',
                              '5' = 'roads',
                              '6' = 'trails'))) %>% 
  
  purrr::set_names('Black bear',
                   'Caribou',
                   'Coyote',
                   'Fisher',
                   'Grey wolf',
                   'Lynx',
                   'Moose',
                   'Red fox',
                   'White-tailed deer')

Merge data into one data frame

coeffs_df_all <- list_rbind(top_models,
                            names_to = 'species') %>% 
  
  mutate(species = as.factor(species),
         
         # add phylopic uuid for each species for plotting 
         # the uuid is extracted using getuuid with the species name as name = ''
         uuid = case_when(species == 'Black bear' ~ get_uuid(name = 'Ursus americanus'),
                          species == 'Caribou' ~ get_uuid(name = 'Rangifer tarandus'),
                          species == 'Coyote' ~ get_uuid(name = 'Canis latrans'),
                          species == 'Fisher' ~ '735066c6-2f3e-4f97-acb1-06f55ae075c9',
                          species == 'Grey wolf' ~ get_uuid(name = 'Canis lupus'),
                          species == 'Lynx' ~ get_uuid(name = 'Lynx lynx'),
                          species == 'Moose' ~ '74eab34a-498c-4614-aece-f02361874f79',
                          species == 'Red fox' ~ '9c977769-bf1e-44d4-82ab-f9f93dce39ca',
                          species == 'White-taield deer' ~ '56f6fdb2-15d0-43b5-b13f-714f2cb0f5d0')) %>% 
  
  # need to remove problematic estimate which is going to skew plot since its so large compared to others
  filter(!c(species == 'Red fox' &
           feature == 'wellsites'))

After plotting the moose image I don’t like it, let’s manually replace it in the data

# I went on the phylopic website and saw there are three images for moose, I like the last one better so we will use it

get_uuid(name = 'Alces alces',
         n = 3)
## [1] "df2d0ad0-adb0-49d7-afe5-edc6cad21064"
## [2] "1a20a65d-1342-4833-a9dd-1611b9fb383c"
## [3] "74eab34a-498c-4614-aece-f02361874f79"
get_uuid(name = 'Odocoileus virginianus',
         n = 3)
## [1] "4584be20-4514-4673-a3e8-97e2a6a10e57"
## [2] "49a5a5db-047e-4e17-849b-9f96a93f0d2b"
## [3] "56f6fdb2-15d0-43b5-b13f-714f2cb0f5d0"
get_uuid(name = 'Pekania pennanti',
         n = 2)
## [1] "5e13bb38-4c9f-4c7a-9c5c-0597fd33e4ab"
## [2] "735066c6-2f3e-4f97-acb1-06f55ae075c9"
get_uuid(name = 'Vulpes',
         n = 2)
## [1] "a1116e25-7b50-4666-bef5-de18b6e2778c"
## [2] "9c977769-bf1e-44d4-82ab-f9f93dce39ca"
# Then I manually copied this uuid and replaces it in the code above

Try plotting all

ggplot(coeffs_df_all, aes(x = feature, 
                          y = estimate, 
                          group = uuid)) +
  
   geom_errorbar(aes(ymin = lower, 
                     ymax = upper, 
                     color = feature),
                width = 0.4,
                linewidth = 0.5,
                position = position_dodge(width = 1.2)) +
  
    # add points for each estimate for each covariate and use position = position_dodge to shift the points so all the species don't plot on top of one another
  geom_phylopic(aes(x = feature, 
                 y = estimate,
                 uuid = uuid),
             position = position_dodge(width = 1.2),
             size = 2)
## Warning: Removed 6 rows containing missing values.
## Warning: `position_dodge()` requires non-overlapping x intervals

ggplot(coeffs, aes(x = feature, y = estimate)) +
  
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
  
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0.4,
                linewidth = 1,
                position = position_dodge(width = 0.5)) +
  
  geom_hline(yintercept = 0, linetype = "dashed")+
  scale_color_manual(values = c("#56B4E9", "#009E73"), name = "Spatial Scale")+
  theme_classic()+
  ggtitle("Moose Response to Anthropogenic Disturbance Features")+
  ylab("Coefficient Estimate \n \u00B1 95% CI")+
  scale_x_discrete(labels =c("Borrowpits", "Harvest\nAreas", "Industrial\nSites", "Pipelines","Roads", "Seismic\nLines", "Trails", "Transmission\nLines"))+
    theme(axis.title.x = element_blank(),
        axis.text.x =  element_text(size = 12),
        axis.title.y = element_text(size = 14),
        legend.title = element_text(size = 12),
        plot.title = element_text(size = 15, hjust = 0.5))

If all else fails can use plot_model function